常微分方程的数值解法实验报告

常微分方程的数值解法
专业班级：信息软件 姓名：吴中原 学号：120108010002 一、实验目的
1、熟悉各种初值问题的算法，编出算法程序；
2、明确各种算法的精度与所选步长有密切关系；通过计算更加了解各种 算法的优越性。

二、实验题目
1、根据初值问题数值算法，分别选择二个初值问题编程计算；
2、试分别取不同步长，考察某节点j x
处数值解的误差变化情况； 3、试用不同算法求解某初值问题，结果有何异常； 4、分析各个算法的优缺点。

三、实验原理与理论基础
（一） 欧拉法算法设计
对常微分方程初始问题
(6-1)
(6-2)
用数值方法求解时，我们总是认为(6-1)、(6-2)的解存在且唯一。

欧拉法是解初值问题的最简单的数值方法。

从(6-2)式由于y (x 0) = y 0已给定，因而可以算出),()('000y x f x y =。

设x 1 = h 充分小，则近似地有：
),()(')
()(00001y x f x y h
x y x y =≈-
(6-3)
记 ,n ,,i x y y i i 10 )(== 从而我们可以取
),(0001y x hf y y ==
作为)(1x y 的近似值。

利用1y 及f (x 1, y 1)又可以算出)(2x y 的近似值：
),(1112y x hf y y +=
一般地，在任意点()h n x n 11+=+处)(x y 的近似值由下式给出
),(1n n n n y x hf y y +=+
(6-4)
这就是欧拉法的计算公式，h 称为步长。

⎪⎩⎪⎨⎧==
)( ),(d d 00y x y y x f x y
（二）四阶龙格-库塔法算法设计：
欧拉公式可以改写为：()
111
,i i i i y y k k hf x y +=+⎧⎪⎨=⎪⎩，它每一步计算(),f x y 的值一次，截
断误差为()2o h 。

改进的欧拉公式可以改写为：()()
()1
1212112
,,i i i i i i y y k k k hf x y k hf x h y k +⎧
=++⎪⎪
=⎨⎪=++⎪⎩
，它每一步要计算(),f x y 的值两次，截断误差为()3o h 。

改进的欧拉方法之所以比欧拉方法具有更高的精度，是因为在每一步它都比欧拉方法多计算了一次(),f x y 的值。

因此，要进一步提高精度，可以考虑在每一步增加计算(),f x y 的次数。

如果考虑在每一步计算(),f x y 的值四次，则可以推得如下公式：
()()()11234121324
31226, 1,221,22,i i
i i i
i i
i i i y y k k k k k hf x y h k hf x y k h k hf x y k k hf x h y k +⎧=++++⎪⎪
⎪=⎪
⎪⎛⎫
=++⎨ ⎪⎝⎭⎪
⎪⎛⎫=++⎪ ⎪⎝⎭⎪
⎪=++⎩ 此公式称为标准四阶龙格-库塔(Runge-Kutta)公式，它的截断误差为()5o h 。

虽然用龙格-库塔方法每一步需要四次调用f ，计算量较改进的欧拉方法大一倍，这里由于龙格-库塔方法的步长增大了一倍，因而两种方法总的计算量相同，但龙格-库塔方法精确度更高。

所以龙格-库塔公式兼顾了精度和计算工作量的较为理想的公式，在实际计算中最为常用。

四、实验内容
（一）问题重述：
科学计算中经常遇到微分方程（组）初值问题，需要利用Euler 法，改进Euler 法，Rung-Kutta 方法求其数值解，诸如以下问题：
（1）()⎪⎩
⎪
⎨
⎧
=-='0 y xy y x y 04 0 < x ≤1
分别取h=0.1,0.2,0.4时数值解。

初值问题的精确解2
54x
e y -+=。

（2）用r=3的Adams 显式和预 - 校式求解
()⎪⎩⎪⎨
⎧-='=-2201y x y y 01≤≤-x
取步长h=0.1，用四阶标准R-K 方法求值。

（3）用改进Euler 法或四阶标准R-K 方法求解
()()()⎪⎩⎪
⎨⎧=-='=-='-='='1033
0012
101213y y y 2y y y y y y
10≤≤x
取步长=h 0.01,计算()()()15.0,10.0,05.0y y y 数值解，参考结果
()()()8613125
.015.0,1493359.015.0,9880787.015.0321≈≈-≈y y y
（4）利用四阶标准R- K 方法求二阶方程初值问题的数值解
（I ）()()⎩
⎨
⎧='==+'-''10,00023y y y y y 0.02h x =≤≤10 (II) ()()⎩⎨
⎧='==+'--''0
0,100)1(1.02y y y y y y 0.1h x =≤≤,10
(III) ()()⎪⎩⎪⎨
⎧='=+='
0,101y y e y y x
0.1h x =≤≤,20 (IV) ()()⎩
⎨
⎧='==+''00,100
sin y y y y 2,400.h x =≤≤
(二)实验代码： 1、欧拉法程序
function y=Euler(a,b,M,y0)
%a=1,b=2,M=10,f=t*y^(1/3),y0=1; h=(b-a)/M;
t=zeros(1,M+1); t=a:h:b;
y=zeros(1,M+1); yy=zeros(1,M+1); y(1)=y0; for k=1:M
y(k+1)=y(k)+h*t(k)*y(k)^(1/3); end
yb=y(M+1);
yy=((t.^2+2)./3).^1.5;
det=yy-y;
plot(t,y,'r-',t,yy,'b:',t,det);
2、改进欧拉法程序
function H=heeuler(a,b,M,ya,f)
%a=0,b=1,M=10,f=t*t+t-y,y0=0;
h=(b-a)/M;
t=zeros(1,M+1);
y=zeros(1,M+1);
p=0;q=0;
t=a:h:b;
y(1)=ya;
for k=1:M
p=feval(f,t(k),y(k));
q=feval(f,t(k+1),y(k)+h*p);
y(k+1)=y(k)+0.5*h*(p+q);
end
yy=t.*t-t+1-exp(-t);
det=yy-y;
plot(t,y,'r-',t,yy,'b:',t,det);
H=[t',y',yy',det']
function f=ff(t,y);
f=t.^2+t-y;
3、四阶龙格-库塔法程序
function H=r_k4(a,b,M,ya,f)
%a=0,b=1,M=10,f=t*t+t-y,y0=0;
h=(b-a)/M;
t=zeros(1,M+1);
t=a:h:b;
y=zeros(1,M+1);
K1=0;K2=0;K3=0;K4=0;
y(1)=ya;
for k=1:M
K1=feval(f,t(k),y(k));
K2=feval(f,t(k)+0.5*h,y(k)+0.5*h*K1);
K3=feval(f,t(k)+0.5*h,y(k)+0.5*h*K2);
K4=feval(f,t(k)+h,y(k)+h*K3);
y(k+1)=y(k)+1/6*(K1+2*K2+2*K3+K4); end
yy=t.*t-t+1-exp(-t);
det=yy-y;
plot(t,y,t,yy,t,det);
H=[t',y',yy',det']
function f=ff(t,y); f=t.^2+t-y;
五、实验结果
1）
()
⎪
⎩
⎪
⎨
⎧
=
-
='
y
xy
y
x
y
4
0 < x≤1
分别取h=0.1,0.2,0.4时数值解。

初值问题的精确解
2
5
4x
e
y-
+
=。

Euler('han',0,0,2,0.1) ans =
0 0
0.1000 0
0.2000 0.0400
0.3000 0.1192
0.4000 0.2356
0.5000 0.3862
0.6000 0.5669
0.7000 0.7729
0.8000 0.9988
0.9000 1.2389
1.0000 1.4874
1.1000 1.7386
1.2000 1.9874
1.3000
2.2289
1.4000
2.4591
1.5000
2.6749
1.6000
2.8736
1.7000 3.0539
1.8000 3.2147
1.9000 3.3561
2.0000
3.4784 Euler('han',0,0,2,0.2) ans =
0 0
0.2000 0
0.4000 0.1600
0.6000 0.4672
0.8000 0.8911
1.0000 1.3886
1.2000 1.9108 1.4000
2.4122 1.6000 2.8568 1.8000
3.2226 2.0000 3.5025 Euler('han',0,0,2,0.4) ans =
0 0 0.4000 0 0.8000 0.6400 1.2000 1.7152 1.6000 2.8119 2.0000 3.5723
2、四阶龙格-库塔法结果
()⎪⎩⎪⎨⎧-='=-2201y x y y 01≤≤-x
取步长h=0.1，用四阶标准R-K 方法求值。

RK4('han',-1,0,1,0.1) ans =
-1.0000 1.0000 -0.9000 0.9909 -0.8000 0.9672 -0.7000 0.9331 -0.6000 0.8921 -0.5000 0.8468 -0.4000 0.7993 -0.3000 0.7515 -0.2000 0.7048 -0.1000 0.6606 0 0.6199
六、实验结果分析与小结
1、步长h 越小则计算精度越高，相对误差越小。

因此，在计算能力允许的范围内，步长越小，得到的结果更精确。

2、改进欧拉法和欧拉法相比较，改进欧拉法的计算精度要更高，相对误差 也较小。

因此在求常微分方程的数值解时，改进欧拉法比欧拉法更精确。