常微分方程初值问题数值解法的比较
数值计算实践—课程设计报告
课题名称
常微分方程初值问题数值解法的比较
完成时间
2013-1-17
姓名
班级
学号
成绩
一. 实验目的及内容
1实验目的：（1）了解常微分方程初值问题的理论背景以及初值问题稳定性、收敛性的研究；
（2）熟练掌握欧拉法、改进欧拉法、龙格-库塔法以及截断误差分析；
输入：[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
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
0.3000 1.3411
0.4000 1.4711
0.5000 1.6082
0.6000 1.7520
0.7000 1.9021
0.8000 2.0580
0.9000 2.2195
1.0000 2.3863
五．计算结果分析:
方法
显示欧拉
简单
精度低
隐式欧拉
稳定性最好
精度低，计算量大
梯形公式
精度提高
计算量大
（4）对右端 用右矩形公式得 ，也叫隐式欧拉法。
误差分析：1.称 为计算 时的局部截断误差；
2.如果数值方法的局部截断误差为 ，那么称这种数值方法的阶数是p，其实p为非负整数。通常情况下，步长h越小，p越高，则局部截断误差越小，计算精
初泰勒展开有
则有 可见欧拉方法是一阶方法，精度不是很高。
2．改进欧拉方法：
考虑一阶常微分方程的初值问题 ，如果存在实数 使得 则称f关于y满足利普希茨条件，L称为利普希茨常数。
对于常微分方程初值问题 ，考虑初值 的扰动是问题的解 发生偏差的情形。若 时 的偏差被控制在有界范围内，则称该初值问题是稳定的，否则该初值问题不稳定的。
特别地，若 时 的偏差收敛于零，则称该初值问题是渐进稳定的。
对于初值问题 稳定性的研究，易知其准确解为 ，假定初值经过扰动后变为 ，对于扰动后的解为 因此带来的扰动误差为 ，因此考虑 时 的值，它取决于 。易知，若 ，则原问题是稳定的；若 ，则原问题是不稳定的；若 ，则原问题是渐进稳定的。
实际遇到的大多数常微分方程初值问题都是稳定的，因此在后面的讨论数值解法时这常常是默认条件。
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);
2.稳定性：
例：考察初值问题 在区间[0,0.5]上的解，分别用欧拉显，隐式格式，改进欧拉格式计算数值解：
节点
欧拉显式
欧拉隐式
改进欧拉法
精确解
0
1.0000
1.0000
1.0000
1.0000
0.1
-2.0000
2.5*10^(-1)
2.5000
4.9787*10^(-2)
0.2
4.0000
6.2500*10^(-2)
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;%根据公式计算
中点公式
精度提高，显式
多一个初值，可能影响精度
1.收敛性：
若某算法对于任意固定的x=xi=x0+i h，当h0 (同时i)时有yiy(xi)，则称该算法是收敛的。
以下讨论的都是单步法（指在计算 时只用到它前一步的信息 ）欧拉法，龙格-库塔法都是单步法的例子；
例：就初值问题 考察欧拉显式格式的收敛性。
二．相关背景知识介绍以及初值问题稳定性的研究：
在科学与工程问题中，常微分方程表述物理量的变化规律，应用非常广泛，比如，天体运动的轨迹，机器人控制，化学反应过程的描述和控制以及电路瞬态过程分析等。这些问题中要求解随时间变化的物理量，即位置函数 表示时间，而微分方程描述了未知函数与它的一阶或高阶导数之间的关系。
6.2500
2.4788*10^(-3)
0.3
-8.0000
1.5625*10^(-2)
1.5625*10^(1)
1.2341*10^(-4)
0.4
1.6000*10^(1)
3.9063*10^(-3)
3.9063*10^(1)
6.1442*10^(-6)
0.5
-3.2000*10^(1)
9.7656*10^(-1)
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.6105 1.7716
2.1436 2.3579 2.5937
1.6165
1.6782
1.7379
龙格-库塔计算结果：
（三阶）输入：Longgekuta3('ff',0,1,1,0.1)
结果：ans =
Longgekuta3('ff',0,1,1,0.1)
ans =
0 1.0000
0.1000 1.1048
0.2000 1.2188
0.3000 1.3411
0.4000 1.4711
end
R=[T' Y'];
四．数值结果：
输入：定义M文件并保存ffx.m
Euler1('ffx',0,1,1,10)
结果：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
梯形公式：对右端 用梯形公式得+显然梯形公式是隐式公式。
改进欧拉公式：先用欧拉公式求的一个初步的近似值 ，成为预测值，预测值 的精度可能达不到要求，在用梯形公式将他校正一次，记为 ，这个结果成为校正值。
预测：
校正：
误差分析：记 为改进欧拉公式在 处的截断误差，
记 因此有
， 表示在 出的局部截断误差。由此得，梯形公式的局部截断误差为 ，因此改进欧拉的截断误差为 ，可见改进欧拉的方法是二阶方法，改进欧拉方法优于欧拉方法。
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;%根据公式计算
向前欧拉公式结果：
输入：Euler2('ffx',0,1,1,10)
结果T =
0 1.0000
0.1000 1.1110
0.2000 1.2344
0.3000 1.3716
0.4000 1.5240
0.5000 1.6933
0.6000 1.8814
0.7000 2.0904
0.8000 2.3227
3．龙格—库塔法:
根据拉格朗日微分中值定理， ，记 得到 ，这样，只要给出一种计算 的算法，就可以得到相应的计算公式。欧拉公式可以写为
改进欧拉公式可以写成 因此推出一般的推出广式
称为p阶龙格-库塔方法，简称p阶R-K方法。
因为 这里的 均为常数。
因为给定的系数不唯一，因此这里的常数有无穷多个解，下面是特殊情况下和一般情况下得结果
（一阶龙格-库塔）当r=1时，这就是欧拉法。
（二阶龙格-库塔）当r=2时， ，这就是改进欧拉法。
三阶和四阶龙格-库塔也只是在一般情况下得结果。
三阶
四阶
其局部截断误差是：
三．程序代码
欧拉法代码如下：
（1）向前欧拉法：
functiony=Euler1(fun,x0,y0,xN,N)
%fun为一阶微分方程，x0，y0为初始条件，xN为取值范围的一个端点，N为区间个数