微分方程数值解上机报告常微分方程初值问题数值解法.docx
《微分方程数值解上机报告常微分方程初值问题数值解法.docx》由会员分享，可在线阅读，更多相关《微分方程数值解上机报告常微分方程初值问题数值解法.docx（9页珍藏版）》请在冰豆网上搜索。
微分方程数值解上机报告常微分方程初值问题数值解法
第一章常微分方程初值问题数值解法
实验一
1、实验题目
试用（a）欧拉格式
（b）中点格式
（c）预报—校正格式
（d）经典四级四阶R-K格式
编程计算方程：
2、程序
#include
#include
#include
constintN=11;
doublefund（doublex,doubley）;
voidEuler（doublea,doubleh,doubley[]）;
voidCenter（doublea,doubleh,doubley[]）;
voidYuJiao（doublea,doubleh,doubley[]）;
voidSjSj（doublea,doubleh,doubley[]）;
voidAdams（doublea,doubleh,doubley[]）;
voidmain（）
{
doublea,b,h,y[N];
inti;
charoption;
cout<<"请输入定义域区间[a,b]：
\n";
cin>>a>>b;
cout<<"请输入初值y[0]：
\n";
cin>>y[0];
h=（b-a）/（N-1）;
cout<<"a:
欧拉法 \nb:
中心法 \nc:
预报-校正格式\n";
cout<<"d:
经典四级四阶R-K格式\n";
cout<<"e:
Adams预报-修正格式 \nf:
退出\n";
label:
cout<<"请选择算法：
";
cin>>option;
switch（option）
{
case'a':
Euler（a,h,y）;break;
case'b':
Center（a,h,y）;break;
case'c':
YuJiao（a,h,y）;break;
case'd':
YuJiao（a,h,y）;break;
case'e':
Adams（a,h,y）;break;
case'f':
exit
（1）;break;
default:
{cout<<"选择有错,重新选择！
\n";gotolabel;}
}
cout<<"计算结果为：
\n xn\t\tyn"<for（i=0;icout<<" "<}
voidEuler（doublea,doubleh,doubley[]）
{
inti;
for（i=1;i{y[i]=y[i-1]+h*fund（a+（i-1）*h,y[i-1]）;}
}
voidCenter（doublea,doubleh,doubley[]）
{
inti;
doublew;
for（i=1;i{w=fund（a+（i-1）*h,y[i-1]）;
y[i]=y[i-1]+h*fund（a+（i-1）*h/2,w*h/2+y[i-1]）;}
}
voidYuJiao（doublea,doubleh,doubley[]）
{
inti;
doublesun;
doubley_Begin[N]={y[0]};
for（i=1;i{y_Begin[i]=y_Begin[i-1]+h*fund（a+（i-1）*h,y_Begin[i-1]）;}
for（i=1;i{
while（fabs（y_Begin[i]-sun）>0.0001）
{
sun=y[i-1]+h/2.0*（fund（a+（i-1）*h,y[i-1]）+fund（a+i*h,y_Begin[i]））;
y_Begin[i]=sun;
}
y[i]=sun;
}
}
voidSjSj（doublea,doubleh,doubley[]）
{
inti;
doublek1,k2,k3,k4;
for（i=1;i{
k1=fund（a+（i-1）*h,y[i-1]）;
k2=fund（a+（i-1）*h+h/2,y[i-1]+h*k1/2）;
k3=fund（a+（i-1）*h+h/2,y[i-1]+h*k2/2）;
k4=fund（a+（i-1）*h+h,y[i-1]+h*k3）;
y[i]=y[i-1]+h/6*（k1+2*k2+2*k3+k4）;
}
}
voidAdams（doublea,doubleh,doubley[]）
{
inti;
doubleme,x[N];
doublek1,k2,k3,k4;
for（i=0;i{
x[i]=a+i*h;
y[i]=y[0];
}
for（i=1;i<4;i++）
{
k1=fund（a+（i-1）*h,y[i-1]）;
k2=fund（a+（i-1）*h+h/2,y[i-1]+h*k1/2）;
k3=fund（a+（i-1）*h+h/2,y[i-1]+h*k2/2）;
k4=fund（a+（i-1）*h+h,y[i-1]+h*k3）;
y[i]=y[i-1]+h/6*（k1+2*k2+2*k3+k4）;
}
for（i=4;i{ me=y[i-1]+h/24*（55*fund（x[i-1],y[i-1]）-59*fund（x[i-2],y[i-2]）+37*fund（x[i-3],y[i-3]）-9*fund（x[i-4],y[i-4]））;
while（fabs（me-y[i]）>0.0001）
{ y[i]=y[i-1]+h/24*（9*fund（x[i],me）+19*fund（x[i-1],y[i-1]）-5*fund（x[i-2],y[i-2]）+fund（x[i-3],y[i-3]））;
me=y[i];
}
}
}
doublefund（doublex,doubley）
{
doubles;
s=y/（2*x）+x*x/2/y;
returns;
}
3、试验结果及分析
（i）运算结果
欧拉法解
中点法解
预报-校正
法解
经典四级四阶R-K法解
1.0
1
1
1
1
1.1
1.1
1.10012
1.1025
1.1025
1.2
1.205
1.20283
1.20997
1.20997
1.3
1.31496
1.3081
1.32235
1.32235
1.4
1.4298
1.4159
1.43954
1.43954
1.5
1.5494
1.52618
1.56141
1.56141
1.6
1.67366
1.6389
1.68785
1.68785
1.7
1.80244
1.75402
1.81873
1.81873
1.8
1.93562
1.87151
1.95393
1.95393
1.9
2.07308
1.99132
2.09334
2.09334
2.0
2.2147
2.11341
2.23683
2.23683
（ii）结果分析
在使用欧拉法数值求解过程中，其计算过程非常简单，即由
可直接计算出
，由
可直接计算出
无需用迭代方法求解任何方程，就可求出近似解
。
通过上面两表的比较，可知随着步长
的减小，近似解
越来越接逼近精确解，即误差越来越小，当步长
充分小时，所得的近似解
能足够地逼近精确解。
通过中点格式
计算出的近似解，格式提高了精度，不过相对于预报-校正格式和经典四级四阶Runge-Kutta格式的高精度，其精度显然还不够高。
在上述两表的比较中，可以看出当步长
取得适当小，用预报—校正格式只需迭代一次就可算出比较好的近似值，得到满足精度要求的解，且比欧拉格式具有更高的精度。
经典四级四阶Runge-Kutta格式是高精度的，并且是显示格式，所以该格式具有算法简单并能得到高精度解的优点。
通过上表可见，四级四阶Runge-Kutta法的确可计算高精度的解。
实验二
1、实验题目
试用Adams预报—修正格式
，
计算下列方程：
3、程序
说明：
本次实验的程序已嵌入到实验一的程序中、下面是本次实验的主要代码。
voidAdams（doublea,doubleh,doubley[]）
{
inti;
doubleme,x[N];
doublek1,k2,k3,k4;
for（i=0;i{
x[i]=a+i*h;
y[i]=y[0];
}
for（i=1;i<4;i++）
{
k1=fund（a+（i-1）*h,y[i-1]）;
k2=fund（a+（i-1）*h+h/2,y[i-1]+h*k1/2）;
k3=fund（a+（i-1）*h+h/2,y[i-1]+h*k2/2）;
k4=fund（a+（i-1）*h+h,y[i-1]+h*k3）;
y[i]=y[i-1]+h/6*（k1+2*k2+2*k3+k4）;
}
for（i=4;i{ me=y[i-1]+h/24*（55*fund（x[i-1],y[i-1]）-59*fund（x[i-2],y[i-2]）+37*fund（x[i-3],y[i-3]）-9*fund（x[i-4],y[i-4]））;
while（fabs（me-y[i]）>0.0001）
{ y[i]=y[i-1]+h/24*（9*fund（x[i],me）+19*fund（x[i-1],y[i-1]）-5*fund（x[i-2],y[i-2]）+fund（x[i-3],y[i-3]））;
me=y[i];
}
}
}
3、试验结果及分析
（i）运行结果
1
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2.0
1
1.1025
1.20996
1.32231
1.43944
1.56125
1.6876
1.81838
1.95346
2.09273
2.23607
（ii）结果分析
升级会员 