常微分方程初值问题数值解法的比较.docx
《常微分方程初值问题数值解法的比较.docx》由会员分享，可在线阅读，更多相关《常微分方程初值问题数值解法的比较.docx（20页珍藏版）》请在冰豆网上搜索。
常微分方程初值问题数值解法的比较
常微分方程初值问题数值解法的比较
数值计算实践—课程设计报告
课题名称
常微分方程初值问题数值解法的比较
完成时间
2013-1-17
姓名
班级
学号
成绩
一.实验目的及内容
1实验目的：
（1）了解常微分方程初值问题的理论背景以及初值问题稳定性、收敛性的研究；
（2）熟练掌握欧拉法、改进欧拉法、龙格-库塔法以及截断误差分析；
（3）比较欧拉法、改进欧拉法及龙格-库塔法，能够选择合适的方法进行问题的研究计算；
2实验内容：
求微分方程（欧拉法求解）
求微分方程（改进欧拉法求解）
求微分方程（龙格-库塔求解）
根据实验的结果进行分析，了解一般方法的的优缺点，稳定性，收敛性以及截断误差的分析，针对相应问题拿出有效方法得出最优的结果。
二．相关背景知识介绍以及初值问题稳定性的研究：
在科学与工程问题中，常微分方程表述物理量的变化规律，应用非常广泛，比如，天体运动的轨迹，机器人控制，化学反应过程的描述和控制以及电路瞬态过程分析等。
这些问题中要求解随时间变化的物理量，即位置函数表示时间，而微分方程描述了未知函数与它的一阶或高阶导数之间的关系。
考虑一阶常微分方程的初值问题，如果存在实数使得则称f关于y满足利普希茨条件，L称为利普希茨常数。
对于常微分方程初值问题，考虑初值的扰动是问题的解发生偏差的情形。
若时的偏差被控制在有界范围内，则称该初值问题是稳定的，否则该初值问题不稳定的。
特别地，若时的偏差收敛于零，则称该初值问题是渐进稳定的。
对于初值问题稳定性的研究，易知其准确解为，假定初值经过扰动后变为，对于扰动后的解为因此带来的扰动误差为，因此考虑时的值，它取决于。
易知，若，则原问题是稳定的；若，则原问题是不稳定的；若，则原问题是渐进稳定的。
实际遇到的大多数常微分方程初值问题都是稳定的，因此在后面的讨论数值解法时这常常是默认条件。
1．欧拉法：
依据：
积分曲线上一点的切线斜率等于函数值。
方法：
推进法，初始点出发，依照方向场在改点的方向推进到
向前欧拉法的得到:
（1）将在处泰勒展开取h的线性部分，得
（2）将初值问题中得导数用向前差商来代替有，因此
（3）将两边同时对x的区间上积分对右端用左矩形公式得，此方法称向前欧拉法，也叫显示欧拉法。
（4）对右端用右矩形公式得，也叫隐式欧拉法。
误差分析：
1.称为计算时的局部截断误差；
2.如果数值方法的局部截断误差为，那么称这种数值方法的阶数是p，其实p为非负整数。
通常情况下，步长h越小，p越高，则局部截断误差越小，计算精
初泰勒展开有
则有可见欧拉方法是一阶方法，精度不是很高。
2．改进欧拉方法：
梯形公式：
对右端用梯形公式得+显然梯形公式是隐式公式。
改进欧拉公式：
先用欧拉公式求的一个初步的近似值，成为预测值，预测值的精度可能达不到要求，在用梯形公式将他校正一次，记为，这个结果成为校正值。
预测：
校正：
误差分析：
记为改进欧拉公式在处的截断误差，
记因此有
，表示在出的局部截断误差。
由此得，梯形公式的局部截断误差为，因此改进欧拉的截断误差为，可见改进欧拉的方法是二阶方法，改进欧拉方法优于欧拉方法。
3．龙格—库塔法:
根据拉格朗日微分中值定理，，记得到，这样，只要给出一种计算的算法，就可以得到相应的计算公式。
欧拉公式可以写为
改进欧拉公式可以写成因此推出一般的推出广式
称为p阶龙格-库塔方法，简称p阶R-K方法。
因为这里的均为常数。
因为给定的系数不唯一，因此这里的常数有无穷多个解，下面是特殊情况下和一般情况下得结果
（一阶龙格-库塔）当r=1时，这就是欧拉法。
（二阶龙格-库塔）当r=2时，，这就是改进欧拉法。
三阶和四阶龙格-库塔也只是在一般情况下得结果。
三阶
四阶
其局部截断误差是：
三．程序代码
欧拉法代码如下：
（1）向前欧拉法：
functiony=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值
end
T=[x',y']
（2）向后欧拉法：
functiony=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;
end
z0=z1;
end
y（n+1）=z1;%%根据向后欧拉公式计算y值
end
T=[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;
fori=1:
n
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;%根据改进欧拉公式计算结果
end
x=x';
y=y';
1．
3.龙格-库塔代码如下：
（三阶龙格-库塔）
functionR=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;%初值赋予
forj=1:
n
k1=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;%根据公式计算
end
R=[T'Y'];
（四阶龙格-库塔）
functionR=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;%初值赋予
forj=1:
n
k1=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;%根据公式计算
end
R=[T'Y'];
四．数值结果：
输入：
定义M文件并保存ffx.m
Euler1（'ffx',0,1,1,10）
结果：
T=
01.0000
0.10001.1000
0.20001.2100
0.30001.3310
0.40001.4641
0.50001.6105
0.60001.7716
0.70001.9487
0.80002.1436
0.90002.3579
1.00002.5937
ans=
Columns1through7
1.00001.10001.21001.33101.46411.61051.7716
Columns8through11
1.94872.14362.35792.5937
向前欧拉公式结果：
输入：
Euler2（'ffx',0,1,1,10）
结果T=
01.0000
0.10001.1110
0.20001.2344
0.30001.3716
0.40001.5240
0.50001.6933
0.60001.8814
0.70002.0904
0.80002.3227
0.90002.5807
1.00002.8674
ans=
Columns1through7
1.00001.11101.23441.37161.52401.69331.8814
Columns8through11
2.09042.32272.58072.8674
改进欧拉公式结果：
输入：
[x,y]=Gaijineuler（'f',0,1,1,0.1）
结果：
x=
0
0.1000
0.2000
0.3000
0.4000
0.5000
0.6000
0.7000
0.8000
0.9000
1.0000
y=
1.0000
1.0959
1.1841
1.2662
1.3434
1.4164
1.4860
1.5525
1.6165
1.6782
1.7379
龙格-库塔计算结果：
（三阶）输入：
Longgekuta3（'ff',0,1,1,0.1）
结果：
ans=
Longgekuta3（'ff',0,1,1,0.1）
ans=
01.0000
0.10001.1048
0.20001.2188
0.30001.3411
0.40001.4711
0.50001.6082
0.60001.7520
0.70001.9021
0.80002.0580
0.90002.2195
1.00002.3863
（四阶）输入：
Longgekuta4（'ff',0,1,1,0.1）
结果：
ans=
01.0000
0.10001.1048
0.20001.2188
0.30001.3411
0.40001.4711
0.50001.6082
0.60001.7520
0.70001.9021
0.80002.0580
0.90002.2195
1.00002.3863
五．计算结果分析:
方法
显示欧拉
简单
精度低
隐式欧拉
稳定性最好
精度低，计算量大
梯形公式
精度提高
计算量大
中点公式
精度提高，显式
多一个初值，可能影响精度
1.收敛性：
若某算法对于任意固定的x=xi=x0+ih，当h0（同时i）时有yiy（xi），则称该算法是收敛的。
以下讨论的都是单步法（指在计算时只用到它前一步的
升级会员 