常微分方程数值解和图形分析资料

第4章常微分方程模型4.3节常微分方程数值解和图形分析4.3.1常微分方程数值解的欧拉方法有很多非线性的常微分方程（组）明明是有解的，可是精确解是非初等函数，不能用基本初等函数经过有限次的四则运算或复合运算来表示，无论是运用在《常微分方程》课程学习过的技巧，还是运用Maple、Mathematica和MATLAB等数学软件，求这些方程（组）的初等函数显式解都是行不通的.4.3.1常微分方程数值解的欧拉方法例如：2dedxyx和ddddxtrxaxyytdybxy研究这些方程（组）的一种替代方法是数值解，就是将微分方程（组）离散化，通过程序计算满足给定的初（边）值条件的解的数值逼近解.微分方程数值解的基石是欧拉方法.4.3.1常微分方程数值解的欧拉方法为了给出微分方程d)d(,fyxxy(4.3.1)的初值问题00dd(,),()yxfxyyxy(4.3.2)在闭区间[a,b]上（其中0ax）的数值逼近解，欧拉（LeonardEuler,1707—1783）引入了迭代公式：11(,)nnnnnnxxhyyhfxy(4.3.3)其中h0，0,1,,()1nbah.（x是不超过x的最大整数）4.3.1常微分方程数值解的欧拉方法欧拉的思想：首先选定水平步长h0，取定实数列{}nx，其中nxanh(0,1,,())nbah.然后用1(,)nnnnyyhfxy作为精确解1()nyx的数值逼近(0,1,,()1)nbah，其中00()yyx.欧拉方法从(xn,yn)到(xn+1,yn+1)的一步(xn,yn)(xn+1,yn)(xn+1,yn+1)hhf(xn,yn)斜率f(xn,yn)(xn+1,y~n+1)图4.174.3.1常微分方程数值解的欧拉方法如图4.17所示，(,)nnfxy是微分方程(4.3.1)式经过点(,)nnxy的积分曲线在点(,)nnxy的斜率，而~1ny是(4.3.1)式满足初始条件()nnyxy的初值问题在1nxx的精确解.欧拉方法从(,)nnxy到11(,)nnxy的一步的局部误差是~1ny与1ny之间的误差，即切线偏离积分曲线的距离.根据泰勒（Taylor）定理可以证明：欧拉方法的局部误差对于h是2阶的，即~211()nnyyOh.4.3.1常微分方程数值解的欧拉方法欧拉方法从00(,)xy到11(,)nnxy的累积误差，就是每一步的局部误差的累积之和，即1ny与初值问题(4.3.2)式的精确解1()nyx之间的误差.定理4.3.1指出欧拉方法的累积误差对于h是1阶的，所以称欧拉方法具有1阶精度.一般的，数值计算方法的精度是由局部误差对于h的阶定义的：如果一个方法的局部误差为1()pOh，就称该方法具有p阶精度.4.3.1常微分方程数值解的欧拉方法定理4.3.1假设初值问题(4.3.2)式在闭区间[a,b]上（其中0ax）有唯一解y(x)，且y(x)在[a,b]上二阶连续可导.如果ny是根据步长h0的欧拉方法在[a,b]上计算得到的对精确值()nyx的逼近值，则存在与h无关的常数C，使得对0,1,,()nbah，都有()nnyyxCh进一步，有0,1,,()0limmax()0nnnbahhyyx.4.3.1常微分方程数值解的欧拉方法欧拉方法(4.3.3)式可以进一步推广为求解常微分方程组初值问题的数值算法.欧拉方法的明显缺点是：（1）理论上，要获得越高的精确性，就需要越小的步长，但计算的时间就会越长；（2）实际上，除了欧拉方法本身的局部误差及累积误差，计算上还存在舍入误差，如果步长太小，舍入误差会累积到不可接受的程度.4.3.1常微分方程数值解的欧拉方法以下是用欧拉方法求解初值问题(4.3.2)式（仅限于一个方程的情形）的MATLAB函数M文件：function[X,Y]=euler(fun,x0,y0,x1,n)%在[x0,x1]计算dy/dx=fun(x,y),y(x0)=y0的数值解h=(x1-x0)./n;X=x0;Y=y0;%将[x0,x1]n等分fork=1:nY(k+1)=Y(k)+h.*fun(X(k),Y(k));X(k+1)=X(k)+h;end4.3.1常微分方程数值解的欧拉方法例4.3.1用欧拉方法计算初值问题2dde,(0)0xyxy在区间[0,1]的数值解.MATLAB脚本：f=@(x,y)exp(-x.^2);[x1,y1]=euler(f,0,0,1,100);[x2,y2]=euler(f,0,0,1,1000);[x3,y3]=euler(f,0,0,1,10000);[x4,y4]=euler(f,0,0,1,100000);[x5,y5]=ode45(f,0:.1:1,0);[0:.1:1;y1(1:10:101);y2(1:100:1001);...y3(1:1000:10001);y4(1:10000:100001);y5.'].'4.3.1常微分方程数值解的欧拉方法表4.4欧拉方法计算初值问题2dde,(0)0xyxyxh=0.01h=0.001h=0.0001h=0.00001ode450000000.10.0997160.0996730.0996680.0996680.0996680.20.197560.197380.197370.197370.197370.30.291660.291280.291240.291240.291240.40.380390.379730.379660.379650.379650.50.462380.461390.461290.461280.461280.60.536660.53530.535170.535160.535150.70.602620.600880.600710.600690.600690.80.660030.657910.657690.657670.657670.90.709010.706520.706270.706240.7062410.749980.747140.746860.746830.746824.3.2常微分方程数值解的MATLAB实现在实践中广泛使用的精度更高的常微分方程数值方法是以德国数学家龙格（Runge）和库塔（Kutta）的姓氏命名的一系列方法，欧拉方法实际上就是1阶龙格-库塔方法.MATLAB实现微分方程数值解的最常用的函数ode45是基于龙格-库塔方法而开发的.MATLAB有7个实现微分方程数值解的函数：ode23、ode45、ode113、ode15s、ode23s、ode23t和ode23tb，它们的语法格式是相同的，其中最常用的是ode45，它是计算常微分方程数值解的首选函数.4.3.2常微分方程数值解的MATLAB实现ode45可以求如下形式的常微分方程（组）初值问题的数值解：00d(,),()dxxxyfyyy(4.3.4)其中T1(),,()nyxyxy，T01,0,0,,nyyy，T1ddd,,dddnyyxxxyT111(,)(,,,),,(,,,)nnnxfxyyfxyyfy积分区间为0[,]fxx或0[,]fxx.4.3.2常微分方程数值解的MATLAB实现ode45的常用语法格式为：（1）[X,Y]=ode45(odefun,xspan,y0)第一输入项odefun是(4.3.4)式的函数f(x,y)的函数句柄，既可以由匿名函数创建，也可以用函数M文件实现，注意odefun的输出必须是列向量.第二输入项xspan是自变量x的要计算数值解的区间，即xspan=[x0,xf]；或者是区间的划分，即xspan=[x0,x1,...,xf]，不要求划分是均匀的，但必须单调增或者单调减.4.3.2常微分方程数值解的MATLAB实现这里的x0和xf即积分区间的0x和fx，不管是0fxx还是0fxx，xspan都必须以0x为第一个元素，以fx为最后一个元素.第三输入项y0即初始条件的0y，虽然0y是列向量，但y0也可以是行向量.4.3.2常微分方程数值解的MATLAB实现第一输出项X是自变量数组，X是列向量.如果xspan=[x0,xf]，则X为ode45迭代计算数值解时产生的自变量x的数列（不一定是等差数列）；如果xspan=[x0,x1,...,xf]，则X=xspan(:)；第二输出项Y是初值问题(4.3.4)式在数组X处的数值解数组，Y的第i列是()iyx的数值解.特别的，如果xspan=[x0,x1,...,xf]，则Y的第i列是()iyx分别在x=x0、x=x1、……、x=xf的数值解.4.3.2常微分方程数值解的MATLAB实现（2）opt=odeset('RelTol',100.*eps,'AbsTol',eps);[X,Y]=ode45(odefun,xspan,y0,opt)第四输入项是将数值解的选项AbsTol（绝对误差精度，默认值1e-6）的值修改成eps（double类型的eps=522≈2.2204e-016），将选项RelTol（相对误差精度，默认值1e-3）修改成100.*eps（系统规定的最小值，如果用户设定的值小于100.*eps，则系统自动设定为100.*eps）.这样做能提高计算精确度，但是需要更长的耗用时间（elapsedtime）.4.3.2常微分方程数值解的MATLAB实现注4.3.1查看MATLAB执行某个（些）语句的耗用时间的方法是：在这个（些）语句的开始前添加命令tic，并在结束后添加命令toc.MATLAB会在命令窗口显示以秒为单位的耗用时间.重复执行同样的语句，会发现耗用时间有一定的随机性，并非每次相同.例4.3.2海上缉私海防某部缉私艇上的雷达发现正东方向c海里处有一艘走私船正以一定速度向正北方向行驶，缉私艇立即以最大速度前往拦截.用雷达进行跟踪时，可保持缉私艇的速度方向始终指向走私船.建立任意时刻缉私艇的位置和缉私艇航线的数学模型，确定缉私艇追上走私船的位置，求出追上的时间.例4.3.2海上缉私1）模型建立和解析解的推导建立直角坐标系，设缉私艇在初始时刻t=0时位于坐标原点，发现走私船位于(c,0)，正以速度a行驶，a的方向始终为y轴正方向，缉私艇立即开始以最大速度b(ba)前往拦截，b的方向始终指向走私船.设缉私艇在任意时刻t位于点P(x,y)，则走私船在同一时刻位于点Q(c,at)，直线PQ与缉私艇航线相切.例4.3.2海上缉私1）模型建立和解析解的推导列式得：22d()d()()xbcxtcxaty(4.3.5)22d()d()()ybatytcxaty(4.3.6)(0)(0)0xy(4.3.7)由(4.3.5)式和(4.3.6)式联立的常微分方程组没有关于x=x(t)和y=y(t)的初等函数显式解.例4.3.2海上缉私1）模型建立和解析解的推导但是可以求得关于y=y(x)的初等函数显式解：将(4.3.6)式除以(4.3.5)式，就可以消去dt以及22()()bcxaty，得d()dycxyatx(4.3.8)例4.3.2海上缉私1）模型建立和解析解的推导为了消去(4.3.8)式的t，注意缉私艇的速度ddbst，其中弧长微分22d(d)(d)sxy，所以2ddd1d1ddddttsyxsxbx(4.3.9)在(4.3.8)式的两边对x求导，并将(4.3.9)式代入，得222dd()1ddyaycxxbx(4.3.10)例4.3.2海上缉私1）模型建立和解析解的推导(4.3.10)式是关于y=y(x)的二阶常微分方程，令ddzyx，可化为关于z=z(x)的一阶常微分方程2d()1dzacxzxb(4.3