《数值分析》教案16.doc

3.4 常微分方程初值问题数值解的MATLAB实现在MATLAB中，有多个求解常微分方程数值解的函数命令，在此介绍几个常用的函数。3.4.1 求常微分方程初值问题数值解的函数ode23和ode45是求解常微分方程数值解最常用的两个函数，命令中的“ode”是英文常微分方程“Ordinary Differential Equation”的缩写，它们都采用龙格-库塔公式进行数值求解，23和45分别表示使用是2/3阶和4/5阶龙格-库塔公式。它们的调用格式基本相同。在此仅以ode23为例来说明函数的用法。函数ode23的调用格式为：x,y=ode23(Fun,Tspan,y0,options)1. 该命令适用于一阶常微分方程组，如遇到高阶常微分方程，必须先把它们变成一阶常微分方程组，即状态方程，方可使用。2. 输入参数“Fun”为定义微分方程组的M-函数文件名，可以在文件名前加写，或用英文格式单引号界定文件名。3. 在编辑调试窗口中编写一阶常微分方程组的M-函数文件时，每个微分方程的格式必须都与一致，即等号左端为待求函数的一阶导数，右端函数的变量严格以“先自变量、后函数”的固定顺序输入，的下标表示微分方程的序数。4. 输入参数“Tspan”规定了常微分方程的自变量取值范围，它以矩阵的形式输入，表示自变量。5. 输入参数“”表示初始条件向量，。微分方程组中的方程个数必须等于初始条件数，这是求常微分方程特解所必须的条件。6. 输入参数“option”表示选项参数，它可由ODESET函数设置，较为复杂。7. 输出参数为微分方程组解函数的列表（x和y都是列矩阵），它包含向量x各节点和与对应向量y的第i个分量值(即第i个方程解)，i表示节点序列数。表3-1 微分方程初值问题的计算函数函数名阶的精确度如何使用ode45中等首先尝试此函数ode23低低精度容差或适度刚性问题ode113从低至高高精度容差ode15s从低至中如果用ode45求解很慢ode23s低低精度容差的线性刚性系统ode23t低问题是适度刚性的ode23tb低用低精度容差求解刚性系统8. 输出参数缺省时输出解函数的曲线，即函数及其各阶导数的曲线。求解微分方程的命令还有ode45，ode113，ode15s，ode23s，ode23t，ode23tb（见表3-1）等。刚性方程组又称为Stiff方程组，其特点是出现解的分量级差别很大，给数值计算带来很大的困难。在化学反应、电子网络和自动控制等领域中经常遇到。3.4.2 ode23与ode45使用方法举例【例3-25】求二阶常微分方程 的通解与特解。解：将二阶常微分方程变换成两个一阶常微分方程组：第一种解法：先建立M-文件。%定义输入,输出变量和函数文件名Function dy=myfun_1(x,y) %明确dy的维数,用微分方程组时不可缺省dy=zeros(2,1); %dy(m)表示y的m阶导数;y(n)表示y的第n列dy(1)=y(2); %与方程组中第二个微分方程对应dy(2)=2*y(2)-2*y(1)+5*exp(2*x)*sin(x); 在命令窗口键入： x y=ode23(myfun_1,0,1,-2;-3)回车得到：x = 0 0.0533 1.0000y = -2.0000 -3.0000 -2.1627 -3.0957 -1.7691 12.8927显示的结果是自变量x和两个待求函数和的对应数据。其中用“”代替了许多输出数据。再键入： x=size(x),y=size(y)x = 15 1y = 15 2表明自变量被分为15个节点，并计算出了对应点上的函数及其一阶导数的取值。第二种解法：使用内联函数inline。在命令窗口键入： myfun_2=inline(y(2);2*y(2)-2*y(1)+5*exp(2*x)*sin(x),x,y); x y=ode23(myfun_2,0 1,-2 -3)回车得到与第一种方法相同的结果。在编辑窗口再键入（不写输出参量）：ode23(myfun_1,0,1,-2;-3)legend(特解函数,一阶导函数) text(0.8,-2,特解函数,FontSize,9)text(0.8,3.5,一阶导函数,FontSize,9)运行得到图3-1所示的微分方程的图示解。图3-1 特解函数及其导数曲线写有输出参量x y时，得出微分方程的数值解，不写输出参量x y时，得出的微分方程的图示解。【例3-26】解微分方程解：用第一种方法解微分方程，建立M-文件：%定义输入,输出变量和函数名function dy=myfun_2(x,y) %明确dy的维数,用微分方程组时不可缺省dy=zeros(3,1); %dy(m)表示y的m阶导数;dy(1)=y(2); %y(n)表示y的第n列dy(2)=y(3); %与方程组中第三个微分方程对应dy(3)=(y(3)-1)2-y(2)-y(1)2; 在命令窗口键入：x y=ode23(myfun_2,0,20,0;1;-1)回车得到： 0 0.0001 19.7666 20.0000 0 1.0000 -1.0000 0.0001 0.9999 -0.9998 0.7802 0.4983 0.2198 0.9015 0.5356 0.0985再键入： x=size(x),y=size(y)x = 105 1y = 105 3表明有105个节点，为函数、一阶导数和二阶导数在每个节点上的取值。第二种解法：使用内联函数inline。在命令窗口键入：myfun_2=inline(y(3);y(2);(y(3)-1)2-y(2)-y(1)2,x,y);x y=ode23(myfun_2,0,20,0;1;-1)回车得到与第一种方法相同的结果。在编辑窗口再键入（不写输出参量）：ode23(myfun_2,0,20,0;1;-1)legend(特解函数,一阶导函数,二阶导函数) text(2,1.6,特解函数,FontSize,9)text(2,0.5,一阶导函数,FontSize,9)text(0.3,-0.7,二阶导函数,FontSize,9)运行得到图3-2所示的微分方程的图示解。写有输出参量x y时，得出微分方程的数值解，不写输出参量x y时，得出的微分方程的图示解。图3-2 特解函数及其一阶、二阶导数曲线【例3-27】求解描述振荡器的经典的Ver der Pol微分方程，。解：编写M-函数文件，设，则：%定义输入,输出变量和函数文件名function dy=myfun_3(t,y) %明确dy的维数,用微分方程组时不可缺省dy=zeros(2,1); %dy(m)表示y的m阶导数;y(n)表示y的第n列dy(1)=y(2); %与方程组中第二个微分方程对应dy(2)=7*(1-y(1)2)*y(2)-y(1); 取，在命令窗口键入：t y=ode45(myfun_3,0,40,1;0)回车得到：t = 0 0.0001 39.9517 40.0000y = 1.0000 0 1.0000 -0.0001 1.7973 -0.1143 1.7918 -0.1149再键入： tn=size(t),yn=size(y)tn = 961 1yn = 961 2表明有961个节点，为函数和一阶导数在每个节点上的取值。下面用两种方法求微分方程的图示解：在命令窗口键入： plot(t,y(:,1),r+,t,y(:,2),g*)legend(特解函数,一阶导函数)text(4,-2.5,特解函数,FontSize,9)text(9,10,一阶导函数,FontSize,9)回车得到图3-3所示的微分方程的图示解。图3-3 与曲线在命令窗口键入： ode45(myfun_3,0,40,1;0)legend(特解函数,一阶导函数) text(4,-2.5,特解函数,FontSize,9)text(9,10,一阶导函数,FontSize,9)回车得到图3-4所示的微分方程的图示解。图3-4 与曲线图3-3与图3-4是用两种不同的方法给出了求微分方程的图示解。【例3-28】解常微分方程Lorenz模型状态方程。若令其初值x0=0;0;1e-10则可以按下面的格式编写出一个MATLAB函数Lorenzq.m来描述系统的动态模型，其内容为function xdot=lorenzq(t,x)xdot=-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3); -x(1)*x(2)+28*x(2)-x(3);在命令窗口键入： ode45(lorenzq,0,100,0;0;1e-10),grid回车得到图3-5所示的微分方程的图示解。图3-5 所示的微分方程的图示解在命令窗口键入： plot(t,x),grid回车得到图3-6状态变量的时间响应曲线。图3-6 状态变量的时间响应曲线在命令窗口键入：plot3(x(:,1),x(:,2),x(:,3),grid,axis(0 40 -20 20 -20 20)回车得到图3-7相空间三维图曲线。图3-7 相空间三维图曲线【例3-29】求微分方程组初值问题选用ode23计算，其相对误差限为，绝对误差限为，分别画出初值条件为=，和解的相平面轨迹图。解：在编辑窗口输入以下命令：ode913=inline(2*x(1)-x(1)*x(2);-x(2)+x(1)*x(2),t,x);option=odeset(reltol,1e-3,abstol,1e-6,outputfcn,odephas2);hold onaxis(0 6 0 6);xlabel(t轴);ylabel(x轴);t,x=ode23(ode913,0 6,1 0.2,option);t,x=ode23(ode913,0 6,1 0.5,option);t,x=ode23(ode913,0 6,1 0.8,option);t,x=ode23(ode913,0 6,1 1.1,option);t,x=ode23(ode913,0 6,1 1.2,option);hold offlegend(初值为1 0.2,初值为1