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

数值分析,数值逼近,实验报告

常微分方程的数值解法

专业班级：信息软件 姓名：吴中原 学号：120108010002

一、实验目的

1、熟悉各种初值问题的算法，编出算法程序；

2、明确各种算法的精度与所选步长有密切关系；通过计算更加了解各种 算法的优越性。

二、实验题目

1、根据初值问题数值算法，分别选择二个初值问题编程计算；

x 2、试分别取不同步长，考察某节点j处数值解的误差变化情况；

3、试用不同算法求解某初值问题，结果有何异常；

4、分析各个算法的优缺点。

三、实验原理与理论基础

（一） 欧拉法算法设计

对常微分方程初始问题

dy f(x,y) dx y(x0) y0 (6-1) (6-2)

用数值方法求解时，我们总是认为(6-1)、(6-2)的解存在且唯一。

欧拉法是解初值问题的最简单的数值方法。从(6-2)式由于y (x0) = y0已给定，因而可以算出y'(x0) f(x0,y0)。

设x1 = h充分小，则近似地有：

y(x1) y(x0) y'(x0) f(x0,y0) h (6-3)

记 yi y(xi) i 0,1, ,n

从而我们可以取

y1 y0 hf(x0,y0)

作为y(x1)的近似值。利用y1及f (x1, y1)又可以算出y(x2)的近似值：

y2 y1 hf(x1,y1)

一般地，在任意点xn 1 n 1 h处y(x)的近似值由下式给出

yn 1 yn hf(xn,yn)

这就是欧拉法的计算公式，h称为步长。 (6-4)

数值分析,数值逼近,实验报告

（二）四阶龙格-库塔法算法设计：

yi 1 yi k1 欧拉公式可以改写为： ，它每一步计算f x,y 的值一次，截 k1 hf xi,yi

断误差为o h2 。 1 y y k1 k2 i i 12 k1 hf xi,yi 改进的欧拉公式可以改写为： ，它每一步要计算

k hf x h,y k ii1 2

f x,y 的值两次，截断误差为o h3 。

改进的欧拉方法之所以比欧拉方法具有更高的精度，是因为在每一步它都比

欧拉方法多计算了一次f x,y 的值。因此，要进一步提高精度，可以考虑在每一步增加计算f x,y 的次数。

如果考虑在每一步计算f x,y 的值四次，则可以推得如下公式：

1y y k1 2k2 2k3 k4 i 1i6 k1 hf xi,yi h1 k hfx ,y k1 2i i22 h1 k3 hf xi ,yi k2 22 k hf x h,y k ii3 4

此公式称为标准四阶龙格-库塔(Runge-Kutta)公式，它的截断误差为o h5 。

虽然用龙格-库塔方法每一步需要四次调用f，计算量较改进的欧拉方法大一倍，这里由于龙格-库塔方法的步长增大了一倍，因而两种方法总的计算量相同，但龙格-库塔方法精确度更高。所以龙格-库塔公式兼顾了精度和计算工作量的较为理想的公式，在实际计算中最为常用。

四、实验内容

（一）问题重述：

科学计算中经常遇到微分方程（组）初值问题，需要利用Euler法，改进Euler 法，Rung-Kutta方法求其数值解，诸如以下问题：

数值分析,数值逼近,实验报告

y 4x

y xy

（1） y 0 0 0 < x 1

分别取h=0.1,0.2,0.4时数值解。

2

初值问题的精确解y 4 5e x。

（2）用r=3的Adams显式和预 - 校式求解

y x2 y2

y 1 0

取步长h=0.1，用四阶标准 1 x 0

R-K 方法求值。

（3）用改进Euler法或四阶标准R-K方法求解

y1 y2 y1 0 1

y y1 y

22 0 0

y 3 y3 y3 0 1 0 x 1

取步长h 0.01,计算y 0.05 , y 0.10 , y 0.15 数值解，

y1 0.15 0.9880787, y2 0.15 0.1493359, y3 0.15 0.8613125

（4）利用四阶标准R- K方法求二阶方程初值问题的数值解

y 3y 2y 0

（I） y 0 0,y 0 1 0 x 1 h 0.02

y 0.1(1 y2)y y 0

(II) y 0 1, y 0 0 0 x 1, h 0.1

y

y

(III) ex 1

y 0 1,y 0 0 0 x 2, h 0.1

y siny 0

(IV) y 0 1,y 0 0 0 x 4, h 0.2

(二)实验代码：

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） 0 < x 1

分别取h=0.1,0.2,0.4时数值解。

xy 4 5e 初值问题的精确解。 24x y xyy y 0 0

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

…… 此处隐藏：2086字，全部文档内容请下载后查看。喜欢就下载吧 ……