数值分析实验报告之常微分方程数值解.docx
《数值分析实验报告之常微分方程数值解.docx》由会员分享，可在线阅读，更多相关《数值分析实验报告之常微分方程数值解.docx（26页珍藏版）》请在冰点文库上搜索。
数值分析实验报告之常微分方程数值解
汐*理^
数学与计算科学学院实验报告
实验项目名称常微分方程数值解
所属课程名称数值方法B
实验
类型验证
实验
日期2013.11.11
班
级
学
号
姓
名
成
绩
-、实验概述：
【实验目的】
1.掌握求解常微分方程的欧拉法；
2.掌握求解常微分方程的预估校正法；
3.掌握求解常微分方程的经典的四阶龙格库塔法；
4.能用C语言或MATLAB将上述三种算法用程序运行出来；
5.将算法实例化，并得出三种算法的相关关系，如收敛性、精度等;
6•附带书中例题的源程序见附录1。
【实验原理】
1•欧拉格式
（1）显式欧拉格式：
yn1=ynhf（Xn,yn）
局部截断误差：
h2h2
y（Xn^Hy^^=—y^-^—y，H（Xn>o（h2）
22
（2）隐式欧拉格式：
yn1二ynhf（Xn1,yn1）
局部截断误差：
h22
y（Xn1）-Yn1：
-『（Xn）=0（h）
2
2•预估校正法
预估：
yn1=ynhf（Xn,yn）
校正：
丄h丄_
yn1=ynAf（Xn,Yn）f（X1,Yn1）]
2
统一格式：
yn^=yn+刁[f（Xn,y.）+f（X^h,『n+hf（X.,y.））]
r
ypin+hf（Xn，yn）,
平均化格式：
丿yc=yn+hf（Xn^yp）,
1=2（yp+yc）.
3•四阶龙格库塔方法的格式（经典格式）
h
yn卅二yn+'Ki+2K2+2K3+K4）,
Ki=f（Xn,Yn）,
hh
K^f（X^-,y^-Ki）,
22
K3=f（Xn十/An+孑K2）,
iK4=f（Xn+h,yn+h&）.
【实验环境】
1.硬件环境：
HP
Microsoft
76481-640-8834005-23929HPCorporation
Intel（R）Core（TM）
I5-2400CPU@3.10GHz
3.09GHz,3.16GB的内存
2.软件环境：
MicrosoftWindowsXP
Professional版本2002ServicePack3
二、实验内容:
【实验方案】
ODE问题:
方案一:
用欧拉法，预估校正法，经典的四阶龙格库塔方法求解下列
例题：
在区间【0,1】上以h=0.1用欧拉法，预估校正法，经典的四阶龙格库塔
法求解微分方程dy/dx=-y+x+1，初值y（0）=1;其精确解为y=x+exp（-x）,且将计
算结果与精确解进行比较，对三个算法的收敛性的进行分析比较。
万案一:
用欧拉法，预估校正法，经典的四阶龙格库塔方法求解初值问题
1
dy/dx=xe"-y，初值y（0）=1;将计算结果与精确解为（x22）e」比较在区
2
间［0,1］上分别取步长h=0.1;0.05时进行计算。
对三个算法的收敛性进行分析比较，
【实验过程】（实验步骤、记录、数据、分析）
注：
以下图形是通过Excel表格处理数据得出，并未通过MATLAB编程序所得!
［业=_y+x+1
1、dx
.y（0）=1
由题可知精确解为：
y=x•/，当x=0时，y（x）=O。
h=0.1
0.8
表1h=0.1时三个方法与精确值的真值表
rH[x.
步长
Euler法
预估校正法
经典四阶库
精确值
0.1
1.010000
1.005000
1.004838
1.249080
0.2
1.029000
1.019025
1.018731
1.055455
0.3[
1.056100
1.041218
1.040818
1.091217
0.4
1.090490
1.070802
1.070320
1.131803
0.5
1.131441
1.107076
1.106531
1.176851
0.6
1.178297
1.149404
1.148812
1.226025
0.7|
1.230467
1.197211
1.196586
1.279016
0.8
1.287421
1.249975
1.249329
1.335536
0.9
1.348678
1.307228
1.306570
1.395322
1.0
1.413811
1.368541
1.367880
1.458127
•Eul日「法
T-预估校正法
Tl经典四阶库
—精确值
图1h=0.1时三个方法走势图
h=0.05（此时将源程序中i的范围进行扩大，即for（i=0;i<20;i++））
表2h=0.05时三个方法与精确值的真值表
rH[x.
步长
Euler法
预估校正法
经典四阶库
精确值
0.05
1.002500
1.001250
1.001229
1.011721
0.10
1.007375
1.004877
1.004837
1.024908
0.15
1.014506
1.010764
1.010708
1.039504
0.20
1.023781
1.018802
1.018731
1.055455
0.25
1.035092
1.028885
1.028801
1.072710
0.30
1.048337
1.040915
1.040818
1.091217
0.35
1.063421
1.054795
1.054688
1.110931
0.40
1.080250
1.070436
1.070320
1.131801
0.45
1.098737
1.087752
1.087628
1.153791
0.50
1.118800
1.106662
1.106531
1.176851
0.55
1.140360
1.127087
1.126950
1.200942
0.60
1.163342
1.148954
1.148812
1.226025
0.65
1.187675
1.172193
1.172046
1.252062
0.70
1.213291
1.196736
1.196585
1.279016
0.75
1.240127
1.222520
1.222367
1.306852
0.80
1.268121
1.249485
1.249329
1.335536
0.85
1.297215
1.277572
1.277415
1.365037
0.90
1.327354
1.306728
1.306570
1.395322
0.95
1.358486
1.336900
1.336741
1.426362
1.00
1.390562
1.368039
1.367880
1.458127
0.8
0.6
0.4
0.2
0
1,6
•Eul^rj去
T-预估枝正法止经典四阶库
―一精确值
00.10.203040.50.6070.80.91
图2h=0.05时三个方法走势图
2、d：
iy
y（0）J
h=0.1
表3h=0.1时三个方法与精确值的真值表
rH[x.
步长
Euler法
预估校正法
经典四阶库
精确值
f0.1
0.900000
0.909625
0.909428
0.929533
0.2
0.819249
0.835927
0.835593
0.872564
10.3
0.754433
0.776081
0.775655
0.826822
[0.4
0.702726
0.727671
0.727189
0.790348
0.5
0.661726
0.688636
0.688127
0.761457
0.6
0.629396
0.657225
0.656711
0.738709
0.7
0.604018
0.631957
0.631453
0.720874
[0.8
0.584147
0.611582
0.611100
0.706908
0.9
0.568575
0.595050
0.594599
0.695927
1.0
0.556297
0.581487
0.581072
0.687191
—•—Euler®
-■■预估校三法
+经典四阶库
精确值
图3h=0.1时三个方法走势图
h=0.05（此时将源程序中i的范围进行扩大，即for（i=0;i<20;i++））
表4h=0.05时三个方法与精确值的真值表
rH[x.
步长
Euler法
预估校正法
经典四阶库
精确值
0.05
0.950000
0.952452
0.952427
0.962924
0.10
0.904904
0.909474
0.909428
0.929533
Qi'「
0.864284
0.870670
0.870606
「0.899511
0.20
0.827741
0.835671
0.835592
0.872564
；0.25
0.794908
0.804137
0.804047
（0.848419
0.30
0.765447
0.775755
0.775655
0.826822
0.35
0.739043
0.750232
0.750125
0.807538
0.40
0.715407
0.727302
0.727189
0.790348
0.45
0.694272
0.706715
0.706599
0.775050
0.50
0.675394
0.688245
0.688126
0.761457
0.55
0.658546
0.671682
0.671561
0.749397
0.60
0.643519
0.656830
0.656710
0.738709
：
0.65:
0.630124
0.643514
0.643395
:
0.729247
0.70
0.618185
0.631570
0.631453
0.720874
0.75
0.607541
0.620848
0.620733
0.713466
[0.80
0.598046
0.611211
0.611100
0.706908
0.85
0.589565
0.602535
0.602426
0.701094
0.90
0.581976
0.594703
0.594599
0.695927
0.95
0.575167
0.587612
0.587512
0.691320
：
1.00
0.569035
0.581167
0.581071
0.687191
・EuIe「法
十预估校正法T-经典四阶库精确值
图4h=0.05时三个方法走势图
【实验结论】（结果）
1.预估校正法的精确度比经典的四阶库法的精确度高，库塔法最低；
2.从表中数据可知三个算法所得数据与精确值相比，可得出以下结论（针对方
案二）：
（1）Euler法所得值偏离精确值最大，因此可知其精度相对来说最差；
（2）经典的四阶库塔法所得值与精确值距离较近，因此可知对于Euler来说,
该法更加有效；
（3）预估校正法的数据时距离精确值最近的，其骤减幅度较小，因此对精度上
的考虑而言，预估校正法应属于最佳解法；
（4）由数据可知，上述两个方程的解的光滑性都比较差，从而导致四阶龙格库塔法的精度低于预估校正法的精度。
3.由图形可知，三个算法所得数据均呈递减趋势，对于他们的收敛性有以下结论：
（1）用上述三法得到的结果大致趋近于0.581，相对于精确值来说，还是存在较大的误差；
（2）就误差最小，应首选预估校正法对问题进行求解，它与经典的四阶库塔法所得值比较接近，因此误差不会相差太大。
【实验小结】（收获体会）
由此次试验，我一方面强化了自己的编程能力，另一方面也对
（1）库塔法，
（2）预估校正法，（3）经典的四阶龙格库塔法有了全新的认识，并能巧妙的将他们运用到数学建模中，解决一些追求高精度的问题。
其次在使用上述三种方法时要充分考虑方程的解的光滑性质，并对照光滑性的
好坏选择相应的求解方法，以达到想要的精度的目的。
三、指导教师评语及成绩：
评语
评语等级
优
良
中
及
格
不及格
1.实验报告按时完成，字迹清楚，文字叙述流畅，逻辑性强
2.实验方案设计合理
3.实验过程（实验步骤详细，记录完整，数据合理，分析透彻）
4实验结论正确.
成绩：
指导教师签名：
批阅日期：
附录1源程序
1.书中例题：
精确值：
#include
#ineludemain（）
{
inti,p;
floatyO,xO,yi,xi,yp,xp,h;xi=0.0;
printf（"请输入步长h:
"）;
scanf（"%f",&h）;for（p=0;pv=10;p++）{
printf（"p=%d",p）;
if（p==O）
{
xp=O.O;
yp=1.0;
}
yp=sqrt（1+2*xp）;
printf（"x%d=%f,y%d=%f\n",p,xp,p,yp）xp+=h;
}
}
Euler格式：
#include
main（）
{
inti,p;
floaty0,x0,yi,xi,yp,xp,h;
xi=0.0;
printf（"请输入步长h:
"）;
scanf（"%f",&h）;
for（i=1;i<=10;i++）
{
p=i-1;
printf（"i=%d",i）;
if（p==0）
{
xp=0.0;
yp=1.0;
}
yi=yp+h*（yp-2*xp/yp）;
printf（"t=%f",xp/yp）;
yp=yi;
xi+=h;
xp=xi;
printf（"x%d=%f,y%d=%f\n",i,xi,i,yi）;}
}
预估校正格式：
#include
main（）
{
inti,t;
floatyi,xi,m,xp,n,yt,xt,h;
xi=O.O;
printf（"请输入步长h:
"）;seanf（"%f",&h）;
printf（"t=OxO=O.OOOOOO,yO=1.OOOOOO\n"）for（i=0;i<10;i++）
{
t=i+1;
printf（"t=%d",t）;
if（i==O）
{
xi=O.O;
yi=1.0;
}
xt=xi+h;
m=yi+h*（yi-2*xi/yi）;
n=yi+h*（m-2*xt/m）;
yt=（m+n）/2.0;
yi=yt;xi+=h;printf（"x%d=%f,y%d=%f\n",t,xt,t,yt）;
}
}
经典的四阶龙格库塔方法的格式：
#include
main（）
{
inti,t;
floatyi,xi,K1,K2,K3,K4,xp,yt,xt,h;
xi=0.0;
printf（"请输入步长h:
"）;scanf（"%f",&h）;
printf（"t=Ox0=0.000000,y0=1.000000\n"）for（i=O;i<1O;i++）
{
t=i+1;
printf（"t=%d",t）;
if（i==0）
{
xi=0.0;yi=1.0;
}
K1=yi-2*xi/yi;
K2=yi+h*K1/2.0-（2*xi+h）/（yi+h*K1/2.0）;
K3=yi+h*K2/2.0-（2*xi+h）/（yi+h*K2/2.0）;
K4=yi+h*K3-2*（xi+h）/（yi+h*K3）;yt=yi+h*（K1+2*K2+2*K3+K4）/6.0;xt=xi+h;
xi+=h;
yi=yt;
printf（"x%d=%f,y%d=%f\n",t,xt,t,yt）;
}
}
2.精确解：
讨e（方案一）
#include
#include
#definee2.1828
main（）
{
inti,p;
floaty0,x0,yi,xi,yp,xp,h;
xi=0.0;
printf（"请输入步长h:
"）;
scanf（"%f",&h）;
for（p=0;p<=10;p++）
{
printf（"p=%d",p）;
if（p==0）
{
xp=0.0;
yp=1.0;
}
yp=xp+pow（e,-xp）;
printf（"x%d=%f,y%d=%f\n",p,xp,p,yp）;xp+=h;
}
}
Euler格式：
#include
main（）
{
inti,p;
floaty0,x0,yi,xi,yp,xp,h;
xi=O.O;
printf（"请输入步长h:
"）;
seanf（"%f",&h）;
for（i=1;i<=11;i++）
{
P=i-1;
printf（"p=%d",p）;
if（p==0）
{
xp=0.0;
yp=1.0;
}
yi=yp+h*（-yp+xp+1）;
yp=yi;
printf（"x%d=%f,y%d=%f\n",p,xi,p,yi）;xi+=h;
xp=xi;
}
}
预估校正格式：
#include
main（）
{
inti,t;
floatyi,xi,m,xp,n,yt,xt,h;
xi=0.0;
printf（"请输入步长h:
"）;
scanf（"%f",&h）;
printf（"t=0x0=0.000000,y0=1.000000\n"）;for（i=0;i<20;i++）
{
t=i+1;
printf（"t=%d",t）;
if（i==0）
{
xi=0.0;
yi=1.0;
}xt=xi+h;
m=yi+h*（-yi+xi+1）;
n=yi+h*（-m+xt+1）;
yt=（m+n）/2.0;
yi=yt;
xi+=h;
printf（"x%d=%f,y%d=%f\n",t,xt,t,yt）;
}
}
经典的四阶龙格库塔方法的格式：
#include
main（）
{
inti,t;
floatyi,xi,K1,K2,K3,K4,yt,xt,h;
xi=0.0;
printf（"请输入步长h:
"）;
scanf（"%f",&h）;
printf（"t=0x0=0.000000,y0=1.000000\n"）;
for（i=0;i<10;i++）
{
t=i+1;
printf（"t=%d",t）;
if（i==0）
{
xi=0.0;
yi=1.0;
}
K1=-yi+xi+1;
K2=-（yi+h*K1/2.0）+xi+h/2.0+1;
K3=-（yi+h*K2/2.0）+xi+h/2.0+1;
K4=-（yi+h*K3）+xi+h+1;yt=yi+h*（K1+2*K2+2*K3+K4）/6.0;xt=xi+h;
xi+=h;
yi=yt;
printf（"x%d=%f,y%d=%f\n",t,xt,t,yt）;
}
}
3.精确解：
1（x22）e"（方案二）
2
#include
#include
#definee2.1828
main（）
{
inti,p;
floatyO,xO,yi,xi,yp,xp,h;
xi=0.0;
printf（"请输入步长h:
"）;
seanf（"%f",&h）;
for（p=0;pv=10;p++）
{
printf（"p=%d",p）;
if（p==0）
{
xp=0.0;
yp=1.0;
}
yp=（xp*xp+2）*pow（e,-xp）/2.0;
printf（"x%d=%f,y%d=%f\n",p,xp,p,yp）;xp+=h;
}
}
Euler格式：
#include
#include
#definee2.1828
main（）
{
inti,p;
floaty0,x0,yi,xi,yp,xp,h;
xi=0.0;
printf（"请输入步长h:
"）;
scanf（"%f",&h）;
printf（"i=0x0=0.000000,y0=1.000000\n"）;for（i=1;i<=10;i++）
{
p=i-1;
printf（"i=%d",i）;
if（p==0）
{
xp=0.0;
yp=1.0;
}
yi=yp+h*（xp*pow（e,-xp）-yp）;
yp=yi;
xi+=h;
xp=xi;
printf（"x%d=%f,y%d=%f\n",i,xi,i,yi）;
}
}
预估校正格式：
#include
#inelude
#definee2.1828
main（）
{
inti,t;
floatyi,xi,m,xp,n,yt,xt,h;
xi=0.0;
printf（"请输入步长h:
"）;
scanf（"%f",&h）;
printf（"i=0x0=0.000000,y0=1.000000\n"）;for（i=0;i<20;i++）
{
t=i+1;
printf（"t=%d",t）;
if（i==0）
{
xi=0.0;
yi=1.0;
}xt=xi+h;
m=yi+h*（xi*pow（e,-xi）-yi）;n=yi+h*（xt*pow（e,-xt）-m）;
yt=（m+n）/2.0;yi=yt;xi+=h;
printf（"x%d=%f,y%d=%f\n",t,xt,t,yt）;
}
}
经典的四阶龙格库塔方法的格式：
#i
升级会员 