Matlab应用12-1.pptVIP

常微分方程 只包含一个自变量的微分方程称常微分方程（Ordinary Differential Equations）。常微分方程没有解析解或是难以求解析解时，应用数值计算求解。 两类：初值问题（Initial Value Problem）和边值问题（Boundary Value Problem） Matlab提供一组配套齐全、结构严谨的指令，包括：微分方程解算（Solver）指令、被解算指令调用的ODE文件格式指令、积分算法参数选项options处理指令，以及输出处理指令等。 解算初值问题的步骤 1、根据具体学科的规律、定理、公式列出微分方程及相应的初始条件 F(y,y’,y’’,…,y(n),t)=0 ; y(0)=y0, y’(0)=y’0,…, y(n-1)(0)=y(n-1)0 y 可以为m*1的列向量。 2、运用变量替换，把高阶方程写成一阶微分方程组，初始条件也做相应变换。 通常令y1=y,y2=y’,…,yn=y(n-1)，原方程改写为： 解算初值问题的步骤 3、编写计算导数的M函数文件(Matlab把它称为ODE文件) 简单方法是把t,Y作为输入宗量，Y’作为输出宗量。也可以包含其它信息。 4、编写好的ODE函数文件和变换后的初值Yo供解算指令调用，得到Y（及其导数）在指定时间区间上的数值解 目前尚没有算法可以有效地解决各种初值问题。为此， Matlab提供了各种解算指令。 与初值常微分方程解算有关的指令表 各种solver解算指令特点 Solver解算指令的使用格式 MATLAB在设计ODE解算软件时，采用相同格式： [t,YY]=solver(‘F’,tspan,Y0) %解算ODE初值问题的最简单调用格式 [t,YY,Te,Ye,Ie]=solver(‘F’,tspan,Y0,options,p1,p2,…) %解算ODE初值问题的最完整调用格式 其中：F是ODE函数文件名（参见例题）；当tspan被赋予二元向量[t0,tf]时，定义求解数值解的时间区间，当tspan被赋予多元向量[t0,t1,…,tf]时， []中元素必须按单调次序排列，求数值解的时间序列；Y0是按照导数阶数递增形式构成的列向量。 输出宗量t是所求数值解的自变量数据列向量(设其长度为N)，而YY是(N×(m×n))矩阵。行中各向量自左向右排列的次序与Y各分量自上而下的次序相同。 options设置算法参数，当输入宗量只有3个或第4宗宗量用空阵赋值时，为缺省设置；p1及以后各输入量是向ODE文件传递参数。 输出宗量中的Te,Ye,Ie是设置“事件(Event)”记录时，才有数据输出。 对微分方程进行解算对应的主程序 global G ME %在主程序中定义全局变量传递参数 1 G=6.672e-11;ME=5.97e24; %输入常量参数G ME vy0=4000;x0=-4.2e7; %输入初值参数vy(0) x(0) t0=0;tf=60*60*24*9;tspan=[t0,tf]; %指定解算微分方程的时间区间 Y0=[x0;0;0;vy0]; %给初值向量赋值 Y(0)=[x(0) 0 0 vy(0)]T [t,YY]=ode45(DYdt,tspan,Y0); %采用4、5阶Runge-Kutta法求解方程 8 X=YY(:,1); %输出YY的第一列，即Y的第一列，对应位移数据x(t) Y=YY(:,2); %输出YY的第二列，即Y的第二列，对应位移数据y(t) plot(X,Y,b,Linewidth,2); hold on axis(image) %x，y坐标轴等长刻度，且坐标框恰包容图形 [XE,YE,ZE] = sphere(15); %产生单位球面数据 RE=0.64e7; %地球半径 XE=RE*XE;YE=RE*YE;ZE=0*ZE; %坐标纸上的地球平面数据 mesh(XE,YE,ZE),hold off % 绘地球示意图 ODE文件和解算指令小结 例题中，Y0和Y的长度是4，所以输出量YY是一个有4列的二维数组，各列依次是[x(t) y(t) vx(t) vy(t)] 编程时注意ODE函数文件[DYdt.m中语句： function Yd=DYdt(t,Y) ] 、解算指令调用格式ode45(‘DYdt’,tspan,Y0);及解算选项(Options)（本例中没有该选项）参数设置3方面的协调。 例题中参数传输通过定