数值计算方法后退欧拉法,龙格库塔法(三阶,四阶方法)

日期： 2011-06-17

     一、    实验目的

1.   学习matlab的使用方法。

2.   掌握常微分方程的几种数值解法：后退欧拉法，龙格库塔法三阶方法，龙格库塔法四阶方法。

3.  比较各方法的数值解及误差，了解各方法的优缺点。

    二、     实验题目

给定的初值问题

  y′=-y+x+2,0≤x≤1

y(0)=-1,

取精确解 y(x)=exp(-x)+x

按   (1) 后退欧拉法，步长h=0.003, h=0.1;

         (2) 龙格库塔法三阶方法，步长h=0.1;

         (3) 龙格库塔法四阶方法，步长h=0.1;

求在节点xk=1+0.1k (k=1,2,3……10)处的数值解及误差比较各方法的优缺点。

    三、   实验原理

1.对于后退欧拉法：

利用Yn+1 = Yn + 1/2*K1 + 1/2*K2①    n = 1, 2, 3……

K1 = hf（Xn, Yn）         ②    

   K2 = hf（Xn + h, Yn + K1） ③    三式可以完成计算

需要将微分方程表达式和精度计算表达式作为两个函数保存在m文件里并在程序中调用：

①微分方程（oulei_wf）

function z=oulei_wf(x,y)

z=-y+x+2

end

②精确解计算(ouleij_q)

function z=ouleij_q(x)

z=exp(-x)+x

end

2.对于龙哥库塔三阶：

利用Yn+1 = Yn + h/6(K1 + 4K2 + K3 )    ①      

              K1 = f（Xn, Yn）                  ②       

              K2 = f（Xn + 1/2*h, Yn +( h/2)*K1）  ③         

    K3= f（Xn + h, Yn - K1+2*h*k2）    ④     四式可以完成计算

3.对于龙哥库塔四阶：

利用Yn+1 = Yn + 1/6(K1 + 2K2 + 2K3 +K4)  ①      

              K1 = hf（Xn, Yn）                  ②       

              K2 = hf（Xn + 1/2*h, Yn + 1/2*K1）   ③       

    K3 = hf（Xn + 1/2*h, Yn + 1/2*K2）   ④       

    K4 = hf（Xn + h, Yn + K3）          ⑤     四式可以完成计算

四、 实验内容

由上述实验原理叙述的后退欧拉法，龙哥库塔三阶，龙哥库塔四阶几种常微分方程数值解法分别对已给定的初值问题进行求解，比较各方法的数值解及误差，了解各方法的优缺点。

五、实验结果

1. 对于后退欧拉法：

①若>>h=0.1;

 y=-1;

 x=1;

for i=1:20

k1=h*oulei_wf(x,y);

k2=h*oulei_wf(x+h,y+k1);

y=y+0.5*k1+0.5*k2

x=x+h;

z=ouleij_q(x)

t=y-z

end

z = 4,z = 3.7000,y = -0.6150

z = 1.4329, z =1.4329, t = -2.0479

z =3.7150, z =3.4435, y = -0.2571

z =1.5012, z =1.5012, t =-1.7583

z =3.4571, z = 3.2114, y =0.0763

z =1.5725, z =1.5725, t =-1.4962

z =3.2237, z = 3.0013, y =0.3876

z =1.66, z =1.66, t =-1.2590

z =3.0124,z =2.8112,y =0.6788

z =1.7231,z =1.7231,t =-1.0444

z =2.8212,z = 2.6391,y =0.9518

z =1.8019,z =1.8019,t =-0.8501

z = 2.82,z =2.4834,y =1.2084

z =1.8827,z =1.8827, t =-0.6743

z =2.4916,z =2.3425, y =1.4501

z =1.9653,z =1.9653, t =-0.5152

2.对于龙哥库塔三阶：

   ①若>> clear all

h=0.1;

 y=-1;

 x=1;

for i=1:25;

k1=oulei_wf(x,y);

k2=oulei_wf(x+h/2,y+k1*(h/2));

k3=oulei_wf(x+h,y-k1*h+k2*(2*h));

y=y+(k1+4*k2+k3)*(h/6);

x=x+h;

z=ouleij_q(x);

t=abs(y-z);

A=[x y z t]

end

z = 4,z = 3.8500,z = 3.7300,z =1.4329

A = 1.1000   -0.6145    1.4329    2.0474

z =3.7145z=3.5788z = 3.4702z =1.5012

A = 1.2000   -0.2562    1.5012    1.7574

z =3.4562z =3.3334z =3.2351z =1.5725

A = 1.3000    0.0776    1.5725    1.4950

z = 3.2224z =3.1113z =3.0224z =1.66

A = 1.4000    0.31    1.66    1.2575

z = 3.0109z = 2.9104z =2.8299z =1.7231

A =1.5000    0.6804    1.7231    1.0427

z =2.8196z =2.7286z =2.6558z =1.8019

A = 1.6000    0.9536    1.8019    0.8483

z =2.z =2.51z =2.4982z =1.8827

A = 1.7000    1.2103    1.8827    0.6724

z =2.47z =2.4152z =2.3556z =1.9653

A =1.8000    1.4521    1.9653    0.5132

z =2.3479z =.2805z =2.2266z =2.0496

A =1.9000    1.6803    2.0496    0.3692

3. 龙哥库塔四阶：

>> clear all

h=0.1;

 y=-1;

 x=1;

for i=1:25;

k1=oulei_wf(x,y);

k2=oulei_wf (x+h/2,y+k1*(h/2));

k3=oulei_wf (x+h/2,y+k2*(h/2));

k4=oulei_wf (x+h,y+h*k3);

y=y+(k1+2*k2+2*k3+k4)*(h/6);

x=x+h;

z=ouleij_q(x);

t=abs(y-z);

A=[x y z t]

end

z =4z =3.8500z =3.8575z =3.7142z =1.4329

A =1.1000   -0.6145    1.4329    2.0474

z =3.7145z =3.5788z =3.5856z =3.4560z =1.5012

A =1.2000   -0.2562    1.5012    1.7574

z =3.4562z =3.3334z =3.3395z =3.2222z =1.5725

A =1.3000    0.0775    1.5725    1.4950

z =3.2225z =3.1113z =3.1169z =3.0108z =1.66

A =1.4000    0.30    1.66    1.2576

z =3.0110z =2.9104z =2.9154z =2.8194z =1.7231

A =1.5000    0.6804    1.7231    1.0427

z =2.8196z =2.7286z =2.7332z =2.63z =1.8019

A =1.6000    0.9536    1.8019    0.8483

z =2.z =2.51z =2.5682z =2.46z =1.8827

A =1.7000    1.2102    1.8827    0.6724

z =2.48z =2.4153z =2.4190z =2.3479z =1.9653

A =1.8000    1.4520    1.9653    0.5133

z =2.3480z =2.2806z =2.2840z =2.2196z =2.0496

A =1.9000    1.6803    2.0496    0.3693

z =2.2197z =2.1587z =2.1618z =2.1035z =2.1353

A =2.0000    1.    2.1353    0.2390

z =2.1036z =2.0485z =2.0512z =1.9985z =2.2225

A =2.1000    2.1014    2.2225    0.1211

z =1.9986z =1.9487z =1.9512z =1.9035z =2.3108

A =2.2000    2.29    2.3108    0.0144

六、实验结果分析

1.对于欧拉法，步长越小，精度越高，而产生的误差越小。总体来说，欧拉法的优点是形式简单，计算方便，缺点是总的运算精度比较低。而且随着x的增大，误差值也越来越大。根据欧拉公式的截断误差计算，欧拉法是一阶方法。

2.对于改进欧拉法，其基本特征与欧拉法相似，也是步长越小，精度越高，误差越小。优点是精度相对欧拉法来说较高，截断误差为O(h^3)，缺点是比欧拉法计算量大。根据欧拉改进法公式的截断误差计算，欧拉改进法是二阶方法。

3.对于龙格—库塔法，优点是精度更高，同样的步长下精度比欧拉法高的多，误差更小 ，截断误差为O(h^5)，故龙格—库塔法是四阶方法。缺点是每步都要计算四次微分值。

4.综上，选取方法时，可综合考虑精度要求和复杂度控制要求等实际需要，从而选择适当的方法求解。