实验二 MATLAB数值计算 常微分方程(组)的求解.doc

实验二 MATLAB数值计算：常微分方程(组)的求解一、实验目的在物理学和工程技术上，很多问题都可以用一个或一组常微分方程来描述，因此要解决相应的实际问题往往需要首先求解对应的微分方程。在大多数情况下这些微分方程通常是非线性的或者是超越方程(比如范德堡方程，波导本征值方程等)，因此往往需要使用计算机数值求解。MATLAB作为一种强大的科学计算语言，其在数值计算和数据的可视化方面具有无以伦比的优势。在解决常微分方程问题上，MATLAB就提供了多种可适用于不同场合(如刚性和非刚性问题)下的求解器(Solver)，例如ode45，ode15s，ode23，ode23s等等。本次实验将以范德堡方程的计算和地球卫星的运行轨道的仿真为例，练习使用MATLAB的常微分方程求解器，以期达到如下几个目的：1. 熟悉常微分方程的求解方法，了解状态方程的概念；2. 能熟练使用dsolve函数解析求解常微分方程；3. 能熟练运用ode45、ode15s求解器分别数值求解非刚性和刚性常微分方程；4. 学习用求解器来绘制相图的方法。二、实验的预备知识1微分方程的概念未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数，称为常微分方程（Ordinary differential equations，简称odes）。n阶常微分方程的一般形式（隐式）为： （1）其中t为自变量。如果未知函数是多元函数，成为偏微分方程。联系一些未知函数的一组微分方程组称为微分方程组。微分方程中出现的未知函数的导数的最高阶解数称为微分方程的阶。若方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，一般表示为若上式中的系数ai (t), i=1,2,n均与t无关，称之为常系数。2常微分方程的解析解有些微分方程可直接通过积分求解，例如一阶常系数常微分方程。有些常微分方程可用一些技巧,如分离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解。线性常微分方程的解满足叠加原理，从而他们的求解可归结为求一个特解和相应齐次微分方程的通解。一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。一阶常微分方程与高阶微分方程可以互化，已知一个n阶常微分方程（显式）： （2）设，可将上式化为n元一阶常微分方程组：（3）注意y1即是要求的未知函数，式(1)、(2)、(3)等价。反过来，在许多情况下，一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根法求解。【(3)式也称状态方程，y1, y2, ,yn称为状态变量，解变量为y 】3微分方程的数值解法除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程主要依靠数值解法。考虑n阶微分方程（2）式的数值求解，它等价于一阶常微分方程组(3)式。将(3)式写成通式：Matlab各种有关的求解器(Solver)都是基于(4)式求解。其中 为n个分量的列向量(Column vector)，也即状态变量； n个分量的列向量，每个元素代表常微分方程组中右边的函数表式 方程(4)式要有确定的解必须给定初值条件(t0为初始时刻)： （7）所谓数值解法，就是在求解区间t0 ,tf上寻求y(t)在一系列离散节点上的近似值。称为步长，通常取为常量h。最简单的数值解法是Euler方法。Euler法的思路极其简单：在节点处用差商近似代替(4)式中的导数利用(4)式导出计算公式（称为Euler格式） （8）即 (9)因此，只要给出初始条件(7),由(8)或(9)式反复迭代，可解出各个时刻tk（k=0,1,m）微分方程的近似解y(tk)。Euler方法只有一阶精度，改进方法有二阶Runge-Kutta法、四阶Runge-Kutta法、五阶Runge-Kutta-Felhberg法和先行多步法等，这些方法可用于解高阶常微分方程（组）初值问题。数值算法的主要缺点是它缺乏物理理解。4解常微分方程的MATLAB命令A. 解析解MATLAB提供了dsolve函数用于计算常微分方程的符号解，其使用格式：S = dsolve (方程1, 方程2,初始条件1，初始条件2 ,自变量) 方程用字符串表示，自变量缺省值为t。导数用D表示，2阶导数用D2表示，以此类推。S返回解析解。在方程组情形，S为一个结构体数组，它的每个域存放方程组的每一个解。例1：求下列微分方程的解析解。 s=dsolve(D2y=sin(2*x)-y,y(0)=0,Dy(0)=1,x); simplify(s) %以最简形式显示sans =-1/3*sin(x)*(-5+2*cos(x) % 结果例2：求微分方程组的解析解。 S=dsolve(Df=f+g,Dg=g-f,Df(0)=1,Dg(0)=1);%默认自变量t simplify(S.f) %S是一个结构体，域f存放解f的结果 simplify(S.g)ans = exp(t)*sin(t) % f(t)ans = exp(t)*cos(t) % g(t)B. 数值解Matlab提供了多种求解器ode23、ode45、ode113、ode15s、ode23s、ode23t 、ode23tb。基本使用格式：tout, Yout = solver (odefun, tspan, y0,options) （solver代表各种求解器，如ode45）变量说明：odefun 一般为用M文件编写的ode函数，须返回一阶常微分方程组(3)式右边的函数值，即f1, f2 , fn。故odefun函数的返回值一般是列向量，其最简单的编写格式为： F = odefun (t, y)t 为标量（自变量），代表迭代点（迭代的当前时刻，如t0, t1等）y 列向量(见5式)，因计算t时刻的函数表达式F需要知道该时刻的y值，见式(8)F 返回值，代表t时刻微分方程组右边的函数值，为列向量（见6式）。有了向量F后，若再给出初始条件，根据(8)或(9)式就可迭代出方程的解。求解器在计算时将会不断地在各个时刻调用odefun的值。此外，MATLAB还专门提供了一个ode函数的编写模板，给出了该函数的详细编写格式，可使用命令: open odefile.m打开，附录2给出了ode函数模板的中文说明。tspan 指定积分迭代的区间t0, tf，t0是迭代的初始时刻。(i). 若 tspan = t0, tf 为2元素的行向量，则求解器的输出矩阵Yout给出求解过程中每个迭代点(即t0, t1, tf)处微分方程组的解，见下面(10)式。(ii). 若tspan = t0, t1, tf为多个元素的行向量（各个时刻点按升序或降序排列），则输出指定时刻点的解。y0 给定初始条件，为n个分量的列向量，见(7)式。options 可选参数，如 RelTol 给出求解方程允许的相对误差，默认10-3AbsTol 给出求解方程允许的绝对误差，默认10-6options参数通过odeset命令设置，odeset 的一般格式如下： options = odeset(属性名1, 属性值1, 属性名2, 属性值2, ) 例如：options = odeset(RelTol,1e-4,AbsTol,1e-4 1e-4 1e-5)。使用help odeset命令可查询求解器的各种参数及其用处。tout 列向量，输出求解过程中各个迭代点的时刻，Yout 矩阵，输出微分方程组(4)式的数值解，每一列对应y的一个分量：Yout的第1列就是要求的未知函数y，第2列是y的一阶导数，第3列是y的二阶导数,提取Y的第1列元素，则可作出yt的图像。ode45是最常用的求解微分方程数值解的命令，对于刚性方程组不宜采用。ode23与ode45类似，只是精度低一些。ode15s用来求解刚性方程组，格式同ode45。可以用help命令查阅有关这些命令的详细信息。表1列出了求解微分方程的各种算法比较：表1：MTALAB提供的各种求解器比较求解器使用算法精度适用ODE类别ode45变步长的四、五阶Runge-Kutta方法适中非刚性的ODEode23二、三阶Runge-Kutta方法低非刚性的ODEode113可变阶的Adams-Bashforth-Moulton PECE算法较高非刚性的ODEode15s可变阶的数值微分方程(NDFs)较低刚性的ODEode23s修正的2阶Rosenbrock方法低刚性的ODEode23t自由插值的梯形法低刚性的ODEode23tb结合梯形法和后向差分公式的Runge-Kutta方法低刚性的ODE说明：1. 所谓刚性方程指的是方程所描写的系统包含多个变化速率差别很大的子系统。如系统：其解 包含快变分量(衰减常数1000)和慢变分量(衰减常数1)，因此上述方程就是刚性的。如果采用非刚性的求解器，那么求解时为了保证误差，积分的步长会很小，从而积分步数很大导致积分时间长的无法进行。2. 在求解ODE方程时，ode45是推荐的首选方法，它是一个中阶的算法。当发现使用ode45解算时效率很低(如计算时间很长)，则应考虑方程可能是刚性的，此时再改用其他的刚性求解器，首选ode15s。3. ode113的求解效率比ode45要好，它适用于对容差有严格要求或一些大型的计算场合。5相图及其绘制以平面相图(二维)为例。考虑一个二阶微分方程，令，可将它转化成状态方程：将方程（11）中的状态变量y1和y2看作是二维平面上的坐标点，这种平面称为状态平面。状态平面上的每一点的坐标(y1，y2)随时间t的变化将在平面上描绘出一条曲线。在给定的初值条件下，方程（11）的解y1(t), y2(t)在平面上绘出的以y1(0), y2(0)为起点的轨迹就称为该状态方程的一个相图（相轨道）。相应地，如果有3个状态变量(y1，y2，y3)则绘出的是三维相图，以此类推。通常从相图可以定性地了解状态方程所描述的系统的工作状态，而不必直接求解微分方程。相图的绘制方法有二： 使用求解指令(如ode45)时，输出变量Yout的每一列代表一个状态变量在不同时刻的值，因此提取Yout每列的值用于绘图，就可绘制出相图。 通过对求解指令的参数进行设置，求解器可以直接输出2维和3维相图，方法是：options = odeset(OutputFcn, odephas3); ode45(odefun,tspan,y0,options);说明： OutputFcn（输出函数）是odeset命令的一个属性名，设置该属性值为odephas3时，则Matlab会自动调用内建函数odephas3用以绘制3维相图；若为odephas2，则自动绘2维相图。三、实验内容1. 求解范德堡电路方程 著名的范德堡电路方程是一个二阶非线性微分方程：令，可将其改写成一阶微分方程组(状态方程)的形式：y1和y2称为状态变量，参量的取值决定了范德堡方程的刚性程度。假定初值条件是：分别两种情况计算范德堡方程的解： =1，计算区间0, 30； =1000，计算区间0, 3000；第一种情况的计算程序如下：global mu; % 定义全局变量，以实现参数在MATLAB的基本工作空间和函数的专用% 空间之间数据的传递mu=1.0; tspan=0,30; % 积分区间y0=0;1; % 初值向量options=odeset(RelTol,1e-03,AbsTol,1e-06,1e-06); %设置参数 % RelTol 相对容许误差，是个标量； % AbsTol 绝对容许误差，标量或矢量。若是矢量则指定y的每个分量(y1,y2)的% 绝对容差；若为标量则指定所有分量的绝对容差t Y=ode45(vandepol,tspan,y0,options); % vandepol是自编ode函数名 % 输出矩阵Y的第一列是范德堡方程的解y(t).%-图形化输出-subplot(2,1,1); % 将图形窗口分成 2x1 的矩阵，第1窗口画y1(t)plot(t,Y(:,1),-o); % Y的第1列是解y(t)xlabel(Time (second),color,r);title(Van de Pol equation for mu=,num2str(mu,5),color,r, FontSize,12); % 加标题并在title中插入变量mu的值,并设置颜色，字号hold on % 添加新的图，将两张图形画在一起plot(t,Y(:,2),-r); % Y的第2列是y2，画y2(t)xlabel(Time (second),color,r);title(Van de Pol equation for mu=,num2str(mu,5),color,r, FontSize,12);legend(y(t),y(t); %加图例：y(t)和y(t), 注意单引号重覆使用代表单引号hold off % 关闭图形添加功能%-下面画相位平面图-subplot(2,1,2); % 第2个窗口画平面相图plot(Y(:,1),Y(:,2); % 画相图xlabel(y1,color,r); ylabel(y2,color,r);title(相平面图,color,r,FontSize,12); % 下面的M函数另存为vandepol.m文件function F=vandepol(t,y)global mu; % 全局变量F=y(2);mu*(1-y(1)2)*y(2)-y(1); % 计算返回值dy/dt，为列向量运行结果如图示：实验任务：就第2种情况编制程序，计算Van de Pol方程的解，并给出结果图和相图。2. 地球卫星运行轨道的仿真(1) 问题的形成轨道上运动的地球卫星，在Newton第二定律，和万有引力定律作用下，有，即这里(x, y)是卫星坐标，是引力常数，是地球的质量。假定卫星以初速度在处进入轨道。令可将(15)式转化为一阶微分方程组：初始条件为：。(2) 实验任务本节练习ode45指令的进阶功能。具体要求如下： 求解方程(16)，绘制相图和卫星轨道图，并标明轨道近地点和远地点，轨道参数（半长轴a，半短轴b，焦距c，偏心率e和轨道周期T）； 在一张图上绘出卫星动能、势能和机械能随时间变化曲线（设卫星质量m=0.5Me）。l 使用sphere()函数绘制3维地球(或卫星)使用格式： X,Y,Z=sphere(N)X,Y,Z = sphere(200); % 产生单位球数据，X,Y,Z均为(200+1)(200+1)的% 矩阵，分别存放球面数据点坐标RE=0.64e7; % 地球半径X=R*X;Y=RE*YE;Z=RE*Z; mesh(X,Y,Z); % 绘制地球l 开启事件判断功能，寻找近地点和远地点任意时刻卫星距离初始点的距离，取它的时间导数du/dt,当du/dt=0时，将给出近地点(u的极小)和远地点(u的极大)坐标(x,y)。当卫星穿越远地点时u的导数值由正向负穿越0；而当卫星穿越近地点时，u的导数值由负向正穿越0，如下图所示。因此可定义一个du/dt的值穿越0的事件，通过穿越0的方式的差别来得到近地点和远地点坐标。在ode函数模板中，case events一项就用于事件的定义和判断，相关的事件子函数events用于编写事件定义以及计算的控制。下面是该事件子函数的程序：function value,isterminal,direction=satellite_event(t,y,y0)% 事件判断子函数。dDSQdt=2*(y(1:2)-y0(1:2)*y(3:4); % dDSQdt 是离初始点的二乘距离 u 的时间导数 du/dt 。 % 当 du/dt=0即dDSQdt=0时触发事件，故以dDSQdt为事件的值value： value = dDSQdt; dDSQdt; % 定义两个穿越0的事件(近地点和远地点事件)：当value值为0时，触发事件终止计算direction = 1; -1; % 第一事件：以渐增方式穿越 0（即从小于0到大于0）(穿越近地点)，对应的(x,y)是% 近地点坐标； % 第二事件：以渐减方式穿越 0(穿越远地点)，对应的(x,y)是远地点坐标； isterminal = 1; 0; % 第一事件发生后，终止计算；而第二事件发生后，继续计算。需要注意的是，要开启事件判断功能，必须设置events参数为on：options=odeset(Events,on);当程序中开启了事件判断功能时，求解器指令的调用格式如下(以ode45为例)：tout, Yout, Te, Ye, Ie = ode45(odefun,tsapn,y0,options,p1,p2,)其中输出变量tout、Yout的含义同前；Te是列向量，存放事件发生的时间；Ye存放事件发生时各个状态变量(这里是y1y4)的值，其每一行代表相应时刻(Te的相应行)各个状态变量的值，因此在本例中Ye的每行头两个元素分别就是近地点和远地点的坐标(x,y)；Ie存放事件的序号，在本例中Ie=1 2 1。输入参数中p1和p2等代表传递的参数，在本例中引力常数G和地球质量Me可以参数的方式传给ode函数(其中计算导数dy/dt时要用)，编程时应尽量避免使用全局变量G和Me，以免引起不必要的麻烦。此外在本例中还有一个参数要传递给ode函数，即初值y0。因为ode函数模板中的事件子函数events(计算近地点和远地点)时需要给出初始位置。ode求解可以有无限多的参数。因此本例的ode指令可以如下书写：options=odeset(Events,on ,OutputFcn,odephas2,);tout, Yout, Te, Ye, Ie = ode45(odefun,tsapn,y0,options,G,Me,y0)其中odefun以用户自编的ode函数名取代，带事件功能的ode函数编写格式可以参考附录2中的模板。附录1：相关的实验结果本次实验的程序可参考附件中的Solve_Stiff_Vandepol.m、vandepol.m和 satellite.m、odefile.m。附录2：MATLAB提供的ode函数编写模板function varargout = odefile(t,y,flag,p1,p2)% 说明：这里给出了ode函数编写的一般格式，函数文件名使用了odefile，用户也可以使用% 其他的文件名。本模板列举的各种情况(case)在具体编写时可依据需要纳入。本ODE函数% 文件自带6个子函数，同样地在编写代码时可依具体情况选用。%t, y分别是自变量(标量)和因变量(列向量)，是最基本的输入宗量(必需)。%flag第三输入宗量，它专供解算指令（如ode45）作调用通知。用户无需直% 接对它赋值，在运行中，解算指令会根据需要向flag发不同的字符串。% p1, p2是被传递的参数。作为示意，仅列出两个。根据需要，传递参数的数目%不受限制。% varargout 是变输出宗量。它由变维的元胞数组构成。每个元胞(即数组元素)中% 可以存放指令产生的任意形式的数据。% 以下是switch-case构成的多向选择控制。应注意：%（A）各情况的情况表达式是MATLAB规定的，不得任意改变。% （B）各情况所执行的任务也是规定的，不得任意改变。%（C）各情况内，变长输出宗量元胞数是规定的，不得任意改变。%（D）各情况内，（等式右边的）函数名称可以改变，但相应子函数名称要一致。%（E）各情况内，（等式右边的）函数中所包含的t , y 是必须的，不管它们是否以显% 式出现相应子函数体中。switch flagcase % 必须用空串。情况为真时，返回 dy/dt = f(t,y)。 varargout1 = f(t,y,p1,p2); % 输出为一个元胞，容纳f子函数的输出dydtcase init % 必须用init，情况为真时，传送计算区间、初值,、设置参数 varargout1:3 = init(p1,p2); % 输出为三个元胞，容纳init子函数的三个输出量case jacobian % 规定jacobian情况：专管计算解析的Jacobian 矩阵 df/dy。 varargout1 = jacobian(t,y,p1,p2);case jpattern % 规定jpattern 情况：专管计算稀疏的数值Jacobian 矩阵 df/dy。 varargout1 = jpattern(t,y,p1,p2);case mass %