常微分方程初值问题数值解法的比较.pdf

1 常微分方程初值问题数值解法的比较数值计算实践 —课程设计报告课题名称常微分方程初值问题数值解法的比较完成时间2013-1-17 姓名班级学号成绩一 . 实验目的及内容1 实验目的：（1）了解常微分方程初值问题的理论背景以及初值问题稳定性、收敛性的研究；（2）熟练掌握欧拉法、改进欧拉法、龙格- 库塔法以及截断误差分析；（3）比较欧拉法、 改进欧拉法及龙格- 库塔法，能够选择合适的方法进行问题的研究计算；2实验内容：求微分方程)10( 1)0(' xyyy（欧拉法求解）求微分方程) 10(1)0(2' xyyxyy（改进欧拉法求解）求微分方程) 10( 1)0() 1log(1' xyxy（龙格 -库塔求解）根据实验的结果进行分析，了解一般方法的的优缺点，稳定性，收敛性以及截断误差的分析，针对相应问题拿出有效方法得出最优的结果二．相关背景知识介绍以及初值问题稳定性的研究：在科学与工程问题中，常微分方程表述物理量的变化规律，应用非常广泛，比如，天体运动的轨迹，机器人控制，化学反应过程的描述和控制以及电路瞬态过程分析等这些问题中要求解随时间变化的物理量，即位置函数tty),(表示时间，而微分方程描述了未知函数与它的一阶或高阶导数之间的关系。

考虑一阶常微分方程的初值问题000')(],,[),,(yxybxxyxfy，如果存在实数,0L使得,,|,|),(|),(|212121RyyyyLyxfyxf则称 f 关于 y 满足 利普希茨条件， L 称为 利普希茨常数2 对于常微分方程初值问题000')(),,(ytyttytfy ， 考虑初值0y的扰动是问题的解)(ty发生偏差的情形若t时)(ty的偏差被控制在有界范围内，则称该初值问题是稳定的 ，否则该初值问题不稳定的 特别地，若t时)(ty的偏差收敛于零，则称该初值问题是渐进稳定的 对于初值问题000')(,ytyttyy 稳定性的研究，易知其准确解为)( 00)(tteyty，假定初值经过扰 动 后 变 为00yy， 对 于 扰 动 后 的 解 为)( 00^0)()(tteyyty因 此 带 来 的 扰 动 误 差 为)( 0^0)(tteyty，因此考虑t时)(ty的值，它取决于易知，若0，则原问题是稳定的；若0，则原问题是不稳定的；若0，则原问题是渐进稳定的实际遇到的大多数常微分方程初值问题都是稳定的，因此在后面的讨论数值解法时这常常是默认条件1． 欧拉法：依据：积分曲线上一点),(yx的切线斜率等于函数值。

方法：推进法，初始点),(000yxp出发，依照方向场在改点的方向推进到nxxx,...,,21向前欧拉法的得到: (1)将)(xy在nx处泰勒展开...! 3)(! 2)()()()()(3' ' '2' ' ' 1hxyhxyhxyxyhxyxynn nnnn取 h 的线性部分，得),(1nnnnyxhfyy（2）将初值问题中得导数用向前差商来代替有),()()()()(111' nnnnnnnnyxfhxyxy xxxyxyy，因此),(1nnnnyxhfyy（ 3 ） 将),(yxfdxdy 两 边 同 时 对x的 区 间],[1nnxx上 积 分111 ))(,()()(即,),(1'nnnnnnxxnnnxxnnxxdxxyxfxyxydxyxfdxy对右端1 ))(,(nnxxndxxyxf用左矩形公式得),(1nnnnyxhfyy，此方法称向前欧拉法，也叫显示欧拉法 4）对右端1 ))(,(nnxxndxxyxf用右矩形公式得),(111nnnnyxhfyy，也叫 隐式欧拉法3 误差分析： 1.称11)(nnyxy为计算1ny时的局部截断误差；2.如果数值方法的局部截断误差为)(1phO，那么称这种数值方法的阶数是p，其实 p 为非负整数。

通常情况下，步长h 越小， p 越高，则局部截断误差越小，计算精nxxy在)(初泰勒展开有)(),(...! 3)(! 2)()()()()(23' ' '2' ' ' 1hOyxhfyhxyhxyhxyxyhxyxynnnnn nnnn则有)()(2 11hOyxynn可见欧拉方法是一阶方法，精度不是很高2． 改进欧拉方法：梯形公式： 对右端1 ))(,(nnxxndxxyxf用梯形公式得+显然梯形公式是隐式公式改进欧拉公式：先用欧拉公式求的一个初步的近似值1ny，成为预测值，预测值1ny的精度可能达不到要求，在用梯形公式将他校正一次，记为1ny，这个结果成为校正值预测：),(1nnnnyxhfyy校正：)],(),([2111nnnnnnyxfyxfhyy误差分析： 记1nT为改进欧拉公式在1nx处的截断误差，)()(11nnnxyxyT记)](,())(,([2)()(1111nnnnnnnxyxfxyxfhxyxyT因此有)])(,())(,([2)()(1111111nnnnnnnnxyxfxyxfhTxyxyT，1nT表示在1nx出的局部 截 断 误 差 由 此 得 ， 梯 形 公 式 的 局 部 截 断 误 差 为)(3 1hOTn， 因 此 改 进 欧 拉 的 截 断 误 差 为)()()(3 11hOxyxyTnnn，可见改进欧拉的方法是二阶方法，改进欧拉方法优于欧拉方法。

3． 龙格—库塔法: 根据拉格朗日微分中值定理，))(,()())(()()(1' 1yhfxyxxyxyxynnnnn，记)),(,(*yhfk得到* 1)()(kxyxynn，这样，只要给出一种计算*k的算法，就可以得到相应的计算公式欧拉公式可以写为 ),(111nnnn yxfKKyy4 改进欧拉公式可以写成),(),()(21121211hKyxfKyxfKKKhyynnnnnn因此推出一般的推出广式),(...),(),()...(1112122122111piipinpnPnnnnPPnnKbhyhaxfKKhbyhaxfKyxfKKCKCKChyy称为 p 阶龙格 -库塔方法，简称p 阶 R-K 方法因为),(),,(,),,(1111ijjijnininnriiinnKhyhxfKyxfKKchyx这里的ijiic,,均为常数因为给定的系数不唯一，因此这里的常数有无穷多个解，下面是特殊情况下和一般情况下得结果（一阶龙格 -库塔）当r=1 时，这就是欧拉法二阶龙格 -库塔）当r=2 时，1,21ijiic，这就是改进欧拉法三阶和四阶龙格-库塔也只是在一般情况下得结果三阶)2,()2,2(),()( 63131213211hKhKyhxfKKhyhxfKyxfKKKKhyynnnnnnnn四阶),()2,2()2,2(),()22(6332212143211hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn其局部截断误差是：111122()()()nnnTy xy xhkk5 三．程序代码欧拉法代码如下：（ 1）向前欧拉法：function y=Euler1(fun,x0,y0,xN,N)%fun 为一阶微分方程，x0 ，y0 为初始条件，xN为取值范围的一个端点，N为区间个数x=zeros(1,N+1);x=zeros(1,N+1);x(1)=x0;y(1)=y0; h=(xN-x0)/N; %h为区间步长for(n=1:N)x(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(fun,x(n),y(n)); % 根据向前欧拉公式计算y 值endT=[x',y']（ 2）向后欧拉法：function y=Euler2(fun,x0,y0,xN,N)% fun为一阶微分方程，x0 ，y0 为初始条件，xN为取值范围的一个端点，N为区间个数 x=zeros(1,N+1);x=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;% h 为区间步长for(n=1:N)x(n+1)=x(n)+h;z0=y(n)+h*feval(fun,x(n),y(n));for(k=1:3)z1=y(n)+h*feval(fun,x(n+1),z0);if(abs(z1-z0)<1e-3)break;endz0=z1;endy(n+1)=z1; %% 根据向后欧拉公式计算y 值endT=[x',y']改进欧拉代码如下：function [x,y]=Gaijineuler(f,x0,y0,xZ,h)%f 为一阶微分方程的函数；x0 ，y0 为初始条件；xZ 为取值范围的一个端点，h为区间步长n=fix((xZ-x0)/h); %计算分点数y(1)=y0;x(1)=x0;for i=1:n6 x(i+1)=x0+i*h;yp=y(i)+h*feval(f,x(i),y(i));yc=y(i)+h*feval(f,x(i+1),yp);y(i+1)=(yp+yc)/2; % 根据改进欧拉公式计算结果endx=x';y=y';1．3. 龙格 -库塔代码如下：（三阶龙格 - 库塔）function R=Longgekuta3(f,a,b,aZ,h)%a，b为端点， h为步长， aZ为初值n=(b-a)/h; T=zeros(1,n+1); %定义向量Y=zeros(1,n+1); T=a:h:b; %计算各个分点Y(1)=aZ; %初值赋予for j=1:nk1=feval(f,T(j),Y(j));k2=feval(f,T(j)+h/2,Y(j)+k1*h/2);k3=feval(f,T(j)+h,Y(j)-h*k1+k2*2*h);Y(j+1)=Y(j)+(k1+4*k2+k3)*h/6; %根据公式计算endR=[T' Y']; （四阶龙格 - 库塔）function R=Longgekuta4(f,a,b,aZ,h)% a ，b为端点， h为步长， aZ 为初值n=(b-a)/h; T=zeros(1,n+1); %定义向量Y=zeros(1,n+1); T=a:h:b; %计算各个分点Y(1)=aZ; %初值赋予for j=1:nk1=feval(f,T(j),Y(j));k2=feval(f,T(j)+h/2,Y(j)+k1*h/2);k3=feval(f,T(j)+h/2,Y(j)+k2*h/2);k4=feval(f,T(j)+h,Y(j)+k3*h);Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)*h/6; %根据公式计算endR=[T' Y'];四．数值结果：输入：定义 M文件并保存 ffx.m Euler1('ffx',0,1,1,10) 7 结果： T = 0 1.0000 0.1000 1.1000 0.2000 1.2100 0.3000 1.3310 0.4000 1.4641 0.5000 1.6105 0.6000 1.7716 0.7000 1.9487 0.8000 2.1436 0.9000 2.3579 1.0000 2.5937 ans = Columns 1 through 7 1.0000 1.1000 1.2100 1.3310 1.4641 1.6。