微分方程数值解上机报告常微分方程初值问题数值解法.docx
《微分方程数值解上机报告常微分方程初值问题数值解法.docx》由会员分享，可在线阅读，更多相关《微分方程数值解上机报告常微分方程初值问题数值解法.docx（29页珍藏版）》请在冰豆网上搜索。
微分方程数值解上机报告常微分方程初值问题数值解法
第一章常微分方程初值问题数值解法
实验一
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<<"计算结果为：
\nxn\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）结果分析
通过实验一的计算结果，可知预报格式已是很好的近似，同时Adams内插法又有较高的精度，故我们用相应的外插公式作为预报格式，然后用内插格式进行迭代，从而形成Adams预报—修正格式，所以Adams预报—修正格式效果显著，能达到较高的精度，尽管只迭代一次却与真解十分相近。
实验三
1、实验题目
分别用四级四阶R-K格式法和Adams预报—较正格式解高阶方程：
2、程序
#include
#include
#include
constm=2;
constn=11;
doubleF1（doublex,doubley1,doubley2）;
doubleF2（doublex,doubley1,doubley2）;
voidmain（）
{
doubleY[m][n],K[m][4];
doubleh,a,b;
inti;
Y[0][0]=0;//初值1
Y[1][0]=1;//初值2
a=0.0;//区间左边界
b=1.0;//区间右边界
h=0.1;//步长
for（i=0;i{
K[0][0]=F1（a+i*h,Y[0][i],Y[1][i]）;
K[1][0]=F2（a+i*h,Y[0][i],Y[1][i]）;
K[0][1]=F1（a+i*h+h/2,Y[0][i]+h*K[0][0]/2,Y[1][i]+h*K[1][0]/2）;
K[1][1]=F2（a+i*h+h/2,Y[0][i]+h*K[0][0]/2,Y[1][i]+h*K[1][0]/2）;
K[0][2]=F1（a+i*h+h/2,Y[0][i]+h*K[0][1]/2,Y[1][i]+h*K[1][1]/2）;
K[1][2]=F2（a+i*h+h/2,Y[0][i]+h*K[0][1]/2,Y[1][i]+h*K[1][1]/2）;
K[0][3]=F1（a+i*h+h,Y[0][i]+h*K[0][2],Y[1][i]+h*K[1][2]）;
K[1][3]=F2（a+i*h+h,Y[0][i]+h*K[0][2],Y[1][i]+h*K[1][2]）;
Y[0][i+1]=Y[0][i]+h/6*（K[0][0]+2*K[0][1]+2*K[0][2]+K[0][3]）;
Y[1][i+1]=Y[1][i]+h/6*（K[1][0]+2*K[1][1]+2*K[1][2]+K[1][3]）;
}
cout<<"用R_K法计算结果为：
\n\n";
cout<<"X\t"<<"y1\t"<<"\ty2\n";
for（i=0;i{
cout<<"x="<cout<cout<<"\n";
}
cout<<"\n";
cout<<"解析精确解为：
\n\n";
cout<<"X\t"<<"y1\t"<<"\ty2\n";
for（i=0;i{
cout<<"x="<cout<cout<<"\n";
}
}
doubleF1（doublex,doubley1,doubley2）
{returny2;}
doubleF2（doublex,doubley1,doubley2）
{return-y1;}
3、试验结果及分析
（i）运行结果
R-K法
精确解
0
0
1
0
1
0.1
0.0998333
0.995004
0.0998334
0.995004
0.2
0.198668
0.980067
0.198669
0.980067
0.3
0.29552
0.955337
0.29552
0.955336
0.4
0.389418
0.921061
0.389418
0.921061
0.5
0.479425
0.877583
0.479426
0.877583
0.6
0.564642
0.825336
0.564642
0.825336
0.7
0.644217
0.764843
0.644218
0.764842
0.8
0.717356
0.696707
0.717356
0.696707
0.9
0.783326
0.621611
0.783327
0.62161
1
0.84147
0.540303
0.841471
0.540302
（ii）结果分析
经典四级四阶Runge-Kutta格式是高精度的，并且是显示格式，所以该格式具有算法简单并能得到高精度解的优点。
通过四级四阶Runge-Kutta格式解高阶方程，由上表可见，四级四阶Runge-Kutta法的确可计算高精度的解。
尽管如此，我们仍不能利用计算机求得它们的精确解，这是因为计算机计算，不论运算器能进行多少位数的运算，总会出现舍入误差以及在计算过程中误差的传递。
有上表的误差分析一栏，可以看出随着计算的进行，误差在不断的累积。
第二章抛物型方程的差分法
实验四
1、实验题目
用古典显式和古典隐式格式计算抛物型方程
满足初始条件
和边界条件
的近似解，
，并与解析解
进行比较（
）。
2、程序
#include
#include
constdoublepi=3.1415926;
constintN=11;
constintM=11;
constdoublet=0.001;
constdoublee=2.718281828459;
doubleUt（doublex）;//初始时刻值
doubleUx1（doubletime）;//左边值
doubleUx2（doubletime）;//右边值
doubleFUN（doublex,doubletime）;
voidmain（）
{
inti,k;
doubleU[M][N];
doubleh,a,b,T1,Tn,r;
cout<<"请输入x所属区间[a,b]\n";
cin>>a>>b;
cout<<"请输入t所属区间（t1,tn）\n";
cin>>T1>>Tn;
h=（b-a）/（N-1）;
r=t/（h*h）;
for（k=0;k{
U[k][0]=Ux1（T1+t*k）;
U[k][N-1]=Ux2（T1+t*k）;
}
for（i=0;iU[0][i]=Ut（a+h*i）;
for（k=1;kfor（i=1;i{U[k][i]=r*U[k-1][i-1]+（1-2*r）*U[k-1][i]+r*U[k-1][i+1];}
cout<<"古典显格式计算结果（t=0.005,h=0.1）：
\n";
for（i=0;i{cout<cout<<"\n";
cout<<"精确计算结果（t=0.005,h=0.1）：
\n";
for（i=0;i{cout<}
doubleUt（doublex）//初始时刻值
{
if（（（int）x）>=1）//为了减小误差当x等于整数时，令x=0；此时sin（x*pi）值不变
x=0;
return（sin（pi*x））;
}
doubleUx1（doubletime）//左边值
{return0;}
doubleUx2（doubletime）//右边值
{return0;}
doubleFUN（doublex,doubletime）
{
doubles1,s2;
if（（（int）x）>=1）
x=0;
s1=-pow（pi,2）*time;
s2=pow（e,s1）*sin（pi*x）;
returns2;
}
3、试验结果及分析
（i）运行结果
节点
古典显式格式解
真解
0
0
0
0.1
0.294186
0.294138
0.2
0.559575
0.559483
0.3
0.770189
0.770063
0.4
0.905411
0.905263
0.5
0.952005
0.95185
0.6
0.905411
0.905263
0.7
0.770189
0.770063
0.8
0.559575
0.559483
0.9
0.294186
0.294138
1
0
0
（ii）结果分析
通过显示差分格式和隐式差分格式求解抛物型方程，由上表的计算结果，可以看出古典隐式格式比古典显式格式具有更高的精度。
实验五
1、实验题目
用Crank-Nicolson差分格式计算抛物型方程
满足初始条件
和边界条件
在
处的解，
。
2、程序
#include
#include
constdoublepi=3.1415926;
constintN=11;
constintM=11;
constdoublet=0.1;
constdoubleh=0.1;
constdoublee=2.71828;
doubleUt（doublex）;//初始时刻值
doubleUx1（doubletime）;//左边值
doubleUx2（doubletime）;//右边值
doubleFUN（doublex,doubletime）;
voidmain（）
{
inti,k;
doubleU[11][11],d[9];
doublea,b,T1,Tn,r;
doubleg[9],w[9];
cout<<"请输入x所属区间[a,b]\n";
cin>>a>>b;
cout<<"请输入t所属区间（t1,tn）\n";
cin>>T1>>Tn;
r=t/（h*h）;
for（k=0;k<11;k++）
{
U[k][0]=Ux1（T1+t*k）;
U[k][10]=Ux2（T1+t*k）;
}
for（i=0;i<11;i++）
U[0][i]=Ut（a+h*i）;
for（k=1;k<11;k++）
{
//计算方程常数项
d[0]=（1-r）*U[k-1][0]+0.5*r*U[k-1][1];
d[8]=（1-r）*U[k-1][9]+0.5*r*U[k-1][8];
for（i=1;i<8;i++）
{
d[i]=（1-r）*U[k-1][i]+0.5*r*（U[k-1][i+1]+U[k-1][i-1]）;
}
g[0]=d[0]/（1+r）;
w[0]=0.5*r/（1+r）;
for（i=1;i<9;i++）
{
g[i]=（d[i]+0.5*r*g[i-1]）/（1+r-0.5*r*w[i-1]）;
w[i]=0.5*r/（1+r-0.5*r*w[i-1]）;
}
U[k][9]=g[8];
for（i=8;i>0;i--）
U[k][i]=g[i-1]+w[i-1]*U[k][i+1];
}
cout<<"古典显格式计算结果（t=0.1,h=0.1）：
\n";
for（i=0;i<11;i++）
{cout<cout<<"\n";
cout<<"精确计算结果（t=0.1,h=0.1）：
\n";
for（i=0;i{cout<cout<<"\n";
cout<<"古典显格式计算结果（t=0.2h=0.1）：
\n";
升级会员 