数值分析实验报告之常微分方程数值解.docx
《数值分析实验报告之常微分方程数值解.docx》由会员分享，可在线阅读，更多相关《数值分析实验报告之常微分方程数值解.docx（24页珍藏版）》请在冰豆网上搜索。
数值分析实验报告之常微分方程数值解
数学与计算科学学院
实验报告
实验项目名称常微分方程数值解
所属课程名称数值方法B
实验类型验证
实验日期2013.11.11
班级
学号
姓名
成绩
一、实验概述：
【实验目的】
1.掌握求解常微分方程的欧拉法；
2.掌握求解常微分方程的预估校正法；
3.掌握求解常微分方程的经典的四阶龙格库塔法；
4.能用C语言或MATLAB将上述三种算法用程序运行出来；
5.将算法实例化，并得出三种算法的相关关系，如收敛性、精度等；
6.附带书中例题的源程序见附录1。
【实验原理】
1．欧拉格式
（1）显式欧拉格式：
局部截断误差：
（2）隐式欧拉格式：
局部截断误差：
2．预估校正法
预估：
校正：
统一格式：
平均化格式：
3．四阶龙格库塔方法的格式（经典格式）
【实验环境】
1.硬件环境：
HP
Microsoft
76481-640-8834005-23929
HPCorporation
Intel（R）Core（TM）
I5-2400CPU@3.10GHz
3.09GHz,3.16GB的内存
2.软件环境：
MicrosoftWindowsXP
Professional
版本2002
ServicePack3
二、实验内容：
【实验方案】
方案一：
用欧拉法，预估校正法，经典的四阶龙格库塔方法求解下列ODE问题：
例题：
在区间【0，1】上以h=0.1用欧拉法，预估校正法,经典的四阶龙格库塔法求解微分方程dy/dx=-y+x+1，初值y（0）=1；其精确解为y=x+exp（-x）,且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较。
方案二：
用欧拉法，预估校正法,经典的四阶龙格库塔方法求解初值问题
dy/dx=
，初值y（0）=1；将计算结果与精确解为
比较在区间[0,1]上分别取步长h=0.1;0.05时进行计算。
对三个算法的收敛性进行分析比较，
【实验过程】（实验步骤、记录、数据、分析）
注：
以下图形是通过Excel表格处理数据得出，并未通过MATLAB编程序所得！
1、
由题可知精确解为：
，当x=0时，y（x）=0。
h=0.1
表1h=0.1时三个方法与精确值的真值表
步长
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
图1h=0.1时三个方法走势图
h=0.05（此时将源程序中i的范围进行扩大，即for（i=0;i<20;i++））
表2h=0.05时三个方法与精确值的真值表
步长
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
图2h=0.05时三个方法走势图
2、
由题可知精确解为：
，当x=0时，y（x）=0。
h=0.1
表3h=0.1时三个方法与精确值的真值表
步长
Euler法
预估校正法
经典四阶库
精确值
0.1
0.900000
0.909625
0.909428
0.929533
0.2
0.819249
0.835927
0.835593
0.872564
0.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
图3h=0.1时三个方法走势图
h=0.05（此时将源程序中i的范围进行扩大，即for（i=0;i<20;i++））
表4h=0.05时三个方法与精确值的真值表
步长
Euler法
预估校正法
经典四阶库
精确值
0.05
0.950000
0.952452
0.952427
0.962924
0.10
0.904904
0.909474
0.909428
0.929533
0.15
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
图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
#include
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=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=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;
}
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=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-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.精确解：
（方案一）
#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=0.0;
printf（"请输入步长h:
"）;
scanf（"%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.精确解：
（方案二）
#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*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
#include
#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）;
}
}
经典的四阶龙格库塔方法的格式：
#include
#include
#definee2.1828
main（）
{
inti,t;
floatyi,xi,K1,K2,K3,K4,yt,xt,h;
xi=0.0;
printf（"请输入步长h:
"）;
scanf（"%f",&h）;
for（i=0;i<10;i++）
{
t=i+1;
printf（"t=%d",t）;
if（i==0）
{
xi=0.0;
yi=1.0;
}
K1=xi*pow（e,-xi）-yi;
K2=（xi+h/2.0）*pow（e,-（xi+h/2.0））-（yi+h*K1/2.0）;
K3=（xi+h/2.0）*pow（e,-（xi+h/2.0））-（yi+h*K2/2.0）;
K4=（xi+h）*pow（e,-（xi+h））-（yi+h*K3）;
yt=yi+h*
升级会员 