数值分析实验报告之常微分方程数值解汇总Word格式.docx
《数值分析实验报告之常微分方程数值解汇总Word格式.docx》由会员分享，可在线阅读，更多相关《数值分析实验报告之常微分方程数值解汇总Word格式.docx（24页珍藏版）》请在冰豆网上搜索。
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时三个方法与精确值的真值表
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
0.25
1.035092
1.028885
1.028801
1.072710
0.30
1.048337
1.040915
0.35
1.063421
1.054795
1.054688
1.110931
0.40
1.080250
1.070436
1.131801
0.45
1.098737
1.087752
1.087628
1.153791
0.50
1.118800
1.106662
0.55
1.140360
1.127087
1.126950
1.200942
0.60
1.163342
1.148954
0.65
1.187675
1.172193
1.172046
1.252062
0.70
1.213291
1.196736
1.196585
0.75
1.240127
1.222520
1.222367
1.306852
0.80
1.268121
1.249485
0.85
1.297215
1.277572
1.277415
1.365037
0.90
1.327354
1.306728
0.95
1.358486
1.336900
1.336741
1.426362
1.00
1.390562
1.368039
图2h=0.05时三个方法走势图
2、
表3h=0.1时三个方法与精确值的真值表
0.900000
0.909625
0.909428
0.929533
0.819249
0.835927
0.835593
0.872564
0.754433
0.776081
0.775655
0.826822
0.702726
0.727671
0.727189
0.790348
0.661726
0.688636
0.688127
0.761457
0.629396
0.657225
0.656711
0.738709
0.604018
0.631957
0.631453
0.720874
0.584147
0.611582
0.611100
0.706908
0.568575
0.595050
0.594599
0.695927
0.556297
0.581487
0.581072
0.687191
图3h=0.1时三个方法走势图
表4h=0.05时三个方法与精确值的真值表
0.950000
0.952452
0.952427
0.962924
0.904904
0.909474
0.864284
0.870670
0.870606
0.899511
0.827741
0.835671
0.835592
0.794908
0.804137
0.804047
0.848419
0.765447
0.775755
0.739043
0.750232
0.750125
0.807538
0.715407
0.727302
0.694272
0.706715
0.706599
0.775050
0.675394
0.688245
0.688126
0.658546
0.671682
0.671561
0.749397
0.643519
0.656830
0.656710
0.630124
0.643514
0.643395
0.729247
0.618185
0.631570
0.607541
0.620848
0.620733
0.713466
0.598046
0.611211
0.589565
0.602535
0.602426
0.701094
0.581976
0.594703
0.575167
0.587612
0.587512
0.691320
0.569035
0.581167
0.581071
图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<
stdio.h>
math.h>
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）;
x%d=%f,y%d=%f\n"
p,xp,p,yp）；
xp+=h;
}
}
Euler格式：
for（i=1;
i++）
p=i-1;
i=%d"
i）;
yi=yp+h*（yp-2*xp/yp）;
t=%f"
xp/yp）;
yp=yi;
xi+=h;
xp=xi;
i,xi,i,yi）;
预估校正格式：
inti,t;
floatyi,xi,m,xp,n,yt,xt,h;
t=0x0=0.000000,y0=1.000000\n"
）；
for（i=0;
10;
t=i+1;
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;
t,xt,t,yt）;
经典的四阶龙格库塔方法的格式：
floatyi,xi,K1,K2,K3,K4,xp,yt,xt,h;
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;
2.精确解：
（方案一）
#definee2.1828
yp=xp+pow（e,-xp）;
p,xp,p,yp）;
xp+=h;
=11;
yi=yp+h*（-yp+xp+1）;
p,xi,p,yi）;
m=yi+h*（-yi+xi+1）;
n=yi+h*（-m+xt+1）;
floatyi,xi,K1,K2,K3,K4,yt,xt,h;
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;
3.精确解：
（方案二）
yp=（xp*xp+2）*pow（e,-xp）/2.0;
i=0x0=0.000000,y0=1.000000\n"
yi=yp+h*（xp*pow（e,-xp）-yp）;
m=yi+h*（xi*pow（e,-xi）-yi）;
n=yi+h*（xt*pow（e,-xt）-m）;
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+
升级会员 