实验五-常微分方程数值解.doc

实验五 常微分方程数值解

———————————————————————————————— 作者：
———————————————————————————————— 日期：

11

个人收集整理 勿做商业用途
实验五 常微分方程数值解
一． 欧拉法

1． 算法说明
对于xi,i=0，1,2,…，n,取步长h为定值时,有h=xi+1-xi，EURLER法的计算公式为：yi+1=yi+h*f(xi,yi）；i=0,1，2,…,n。

2． 程序中主要符号说明
a为x的下界,
b为x的上界，
h为步长，
n为循环次数（即为x的数值点数减一），
x0、y0为循环初始值，
x1、y1为输出值。

3． 计算流程框图
开始
读入数据x0，y0，h，n
FOR I=1 TO N
X0+H=>X1,Y0+H*F(X0,Y0)=Y1
输出X1，Y1
X1=>X0
Y1=>Y0
NEXT
结束

4． 程序清单
syms a b h n x0 y0 x1 y1；
a=0;
b=0.1；
n=5;
x0=0;
y0=1;
h=（b—a)/n；
for j=1:n
x1=x0+h;
y1=y0-h＊0。9*y0/(1+2＊x0);
x1
y1
x0=x1;
y0=y1；
end

5． 算例及输入数据说明
算例：求解初值问题

当x=0，0.02,0.04，…，0.10时的数值解。

输入数据说明:
a=0； a为x 的下界
b=0。1； b为x的上界
n=5; n为循环次数,即为x的数值点数减一
x0=0； x（0）的值
y0=1； y（0)的值

6． 程序运行结果及结果分析
运行结果：
x1 =

0.0200
y1 =

0.9820
x1 =

0.0400
y1 =

0.9650
x1 =

0。0600
y1 =

0。9489
x1 =

0.0800
y1 =

0。9337
x1 =

0。1000
y1 =

0.9192
结果分析:欧拉法计算简单，但计算效率并不高，计算精度很低，局部截断误差较大.

二． 改进欧拉法

1． 算法说明
对于xi,i=0,1,2,…，n，取步长h为定值时，有h=xi+1-xi，EURLER法的计算公式为:　yp=yi+h＊f(xi，yi）
　　 yc=yi+h*f(xi+1，yp）
　　 yi+1=（yp+yc)/2；i=0,1,2，…，n。

2． 程序中主要符号说明
a为x的下界，
b为x的上界，
h为步长，
n为循环次数（即为x的数值点数减一)，
x0、y0为循环初始值，
yp、yc为运算中间值，
x1、y1为输出值.

3． 计算流程框图
开始
读入数据x0，y0，h，n
FOR I=1 TO N
X0+H=>X1 Y0+H*F(X0,Y0)=>YP
Y0+H*F(X0,Y0)=>YC (YP+YC)/2=>Y1

输出X1，Y1
X1=>X0
Y1=>Y0
NEXT
结束

4． 程序清单
syms a b h n x0 y0 yp yc x1 y1；
a=0;
b=0。1;
n=5；
x0=0;
y0=1;
yp=1;
yc=1；
h=(b—a)/n;
for j=1：n
x1=x0+h;
yp=y0—h＊0。9*y0/(1+2*x0)；
yc=y0—h＊0.9＊yp/（1+2*x1);
y1=（yp+yc)/2;
x1
y1
x0=x1;
y0=y1;
end

5． 算例及输入数据说明
算例：求解初值问题

当x=0，0。02,0.04，…，0。10时的数值解。

输入数据说明:
a=0； a为x的下界
b=0.1； b为x的上界
n=5； n为循环次数，即为x的数值点数减一
x0=0; x(0)的值
y0=1; y（0）的值
yp=1; 由y(0)的值决定
yc=1; 由y（0）的值决定

6． 程序运行结果及结果分析
x1 =

0.0200
y1 =

0.9825
x1 =

0。0400
y1 =

0。9660
x1 =

0.0600
y1 =

0.9503
x1 =

0。0800
y1 =

0。9354
x1 =

0。1000
y1 =

0。9212

结果分析:
计算过程比欧拉法较复杂，但改进欧拉法先用欧拉法求出预报值,再利用公式求出校正值，局部截断误差比欧拉法低了一阶，较大程度地提高了计算精度.

三． 龙格库塔法

1． 算法说明
对于xi，i=0，1，2，…，n，取步长h为定值时，有h=xi+1-xi,EURLER法的计算公式为：yi+1＝yi＋h＊( K1+ 2＊K2 +2*K3+ K4）/6
　　K1=f(xi,yi)
　　K2=f(xi+h/2,yi+h*K1/2）
　　K3=f（xi+h/2,yi+h*K2/2）
　　K4=f(xi+h,yi+h*K3); i=0,1，2，…，n。

2． 程序中主要符号说明
h为步长,
n为循环次数（即为x的数值点数减一)，
x0、y0为循环初始值,
k1，k2，k3,k4为运算中间值，
x1、y1为输出值。

3． 计算流程框图
开始
读入数据x0，y0，h，n
FOR I=1 TO N
X0+H=>X1 F(X0,Y0)=>K1
F(X0+H/2,Y0+K1*H/2)=>K2
F(X0+H/2,Y0+K2*H/2)=>K3
F(X0+H,Y0+K3*H)=>K4
Y0+(K1+2K2+2K3+K4)=>Y1

输出X1，Y1
X1=>X0
Y1=>Y0
NEXT
结束

4． 程序清单
syms h n x0 y0 x1 y1 k1 k2 k3 k4;
x0=0;
y0=1；
h=0.2;
n=5;
for j=1：n
x1=x0+h;
k1=x0+y0；
k2=x0+h/2+y0+k1*h/2;
k3=x0+h/2+y0+k2*h/2;
k4=x0+h+y0+k3＊h；
y1=y0+h＊(k1+2＊k2+2*k3+k4)/6;
x1
y1
x0=x1；
y0=y1;
end
5． 输入数据说明
算例：求解初值问题

输入数据说明：
x0=0； x(0）的值
y0=1； y（0）的值
h=0.2; 步长为0.2
n=5; n为循环次数，即为x的数值点数减一

6． 程序运行结果及结果分析
x1 =

0.2000
y1 =

1.2428
x1 =

0.4000
y1 =

1.5836
x1 =

0。6000
y1 =

2.0442
x1 =

0.8000
y1 =

2.6510
x1 =

1
y1 =

3.4365

数据分析:
龙格-库塔法是应用广泛的高精度单步算法。此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂。