微分方程第一章习题5实验报告讲解.docx
《微分方程第一章习题5实验报告讲解.docx》由会员分享，可在线阅读，更多相关《微分方程第一章习题5实验报告讲解.docx（22页珍藏版）》请在冰豆网上搜索。
微分方程第一章习题5实验报告讲解
第一章习题5实验报告
一、实验目的
通过用Euler法、中点格式、预报－校正格式、Adams预报—修正格式等算法求解一阶常微分方程初值问题
的数值解，经典四级四阶R-K法是高精度单步法，四阶Adams预估－校正法是线性多步法，这两个算法较实验一算法精度要高。
通过本次实验，掌握R-K法的编程要领，掌握多步法初值的准备方法，深刻领会微分方程数值解的实质，体会单步法和线性多步法各自的优缺点，熟练掌握各算法的计算机实现过程，并能从理论及实验结果分析各种算法的优缺点。
二、实验内容
用Euler法、中点格式、预报－校正格式、经典四级四阶R-K格式、Adams预报—修正格式算法求解一阶常微分方程初值问题
的数值解，其精确解为
。
三、算法
1Euler法
其中
，
为步长。
2 中点格式
3 预报－校正格式
4四级四阶格式：
5Adams预报—修正格式
其中
为步长。
四程序
一、Euler法程序如下：
#include
#defineLEN100
intEuler（inta,intb,doubleh）
{
doubley[LEN]={0};
doublex[LEN]={0};
inti=0;
doublenum=0;
num=（（b-a）/h）+1;//节点个数
//循环计算出各个节点的坐标
for（i=1;i{
x[0]=a;
x[i]=a+i*h;
}
//带公式计算各个节点的y的值
for（i=0;i{
y[0]=1;
y[i+1]=（y[i]/（2*x[i]））+（x[i]*x[i]）/（2*y[i]）;
}
//输出x和y
printf（"\n*************************Euler法求解方程*******************************\n"）;
printf（"\nx[n]y[n]\n\n"）;
for（i=0;i{
printf（"%.4f%.4f\n",x[i],y[i]）;
}
printf（"\n*************************Euler法求解方程********************************\n\n"）;
return0;
}
}
二、中点格式程序如下：
intMid（inta,intb,doubleh）
{
doubley[LEN]={0};
doublex[LEN]={0};
inti=0;
doublenum=0;
num=（（b-a）/h）+1;//节点个数
//循环计算出各个节点的坐标
for（i=0;i{
x[0]=a;
x[i]=a+i*h;
}
//带公式计算各个节点的y的值
for（i=0;i{
y[0]=1;
y[i+1]=y[i]+h*f（（x[i]+h/2）,（y[i]+（h/2）*（f（x[i],y[i]））））;
}
//输出x和y
printf（"\n*************************中点格式求解方程*******************************\n"）;
printf（"x[n]y[n]\n\n"）;
for（i=0;i{
printf（"%.4f%.4f\n",x[i],y[i]）;
}
printf（"\n*************************中点格式求解方程*******************************\n\n"）;
return0;
}
三、预报——校正格式程序如下：
intPreCorr（inta,intb,doubleh）
{
doubley1[LEN]={0};
doubley[LEN]={0};
doublex[LEN]={0};
inti=0;
doublenum=0;
num=（（b-a）/h）+1;//节点个数
//循环计算出各个节点的坐标
for（i=0;i{
x[0]=a;
x[i]=a+i*h;
}
//带公式计算各个节点的y的值
for（i=0;i{
y[0]=1;//初值为1
y1[i+1]=y[i]+h*f（x[i],y[i]）;
y[i+1]=y[i]+（h/2）*（f（x[i],y[i]）+f（x[i+1],y1[i+1]））;
}
//输出x和y
printf（"\n***********************预报—校正格式求解方程***************************\n"）;
printf（"x[n]y[n]\n\n"）;
for（i=0;i{
printf（"%.4f%.4f\n",x[i],y[i]）;
}
printf（"\n***********************预报—校正格式求解方程***************************\n\n"）;
return0;
}
四、经典四级四阶R-K程序如下：
intRungeKutta（inta,intb,doubleh）
{
doubley1[LEN]={0};
doubley[LEN]={0};
doublex[LEN]={0};
inti=0;
doublenum=0;
num=（（b-a）/h）+1;//节点个数
//循环计算出各个节点的坐标
for（i=0;i{
x[0]=a;
x[i]=a+i*h;
}
//带公式计算各个节点的y的值
for（i=0;i{
doubleK1=0;
doubleK2=0;
doubleK3=0;
doubleK4=0;
y[0]=1;//初值为1
K1=f（x[i],y[i]）;
K2=f（（x[i]+h/2）,（y[i]+h*K1/2））;
K3=f（（x[i]+h/2）,（y[i]+h*K2/2））;
K4=f（（x[i]+h）,（y[i]+h*K3））;
y[i+1]=y[i]+（h/6）*（K1+2*K2+2*K3+K4）;
}
//输出x和y
printf（"\n********************经典四级四阶R-K格式求解方程*************************\n"）;
printf（"x[n]y[n]\n\n"）;
for（i=0;i{
printf（"%.4f%.4f\n",x[i],y[i]）;
}
printf（"\n********************经典四级四阶R-K格式求解方程*************************\n\n"）;
return0;
}
六、Adams程序如下：
intAdams（inta,intb,doubleh）
{
doubley1[LEN]={0};
doubley2[LEN]={0};
doubley[LEN]={0};
doublex[LEN]={0};
inti=0;
doublenum=0;
doubletmp=0;
num=（（b-a）/h）+1;//节点个数
//循环计算出各个节点的坐标
for（i=0;i{
x[0]=a;
x[i]=a+i*h;
}
for（i=0;i{
y2[i]=sqrt（x[i]/2+（x[i]*x[i]）/2）;
}
y2[0]=1;
//带公式计算各个节点的y的值
for（i=0;i<4;i++）
{
doubleK1=0;
doubleK2=0;
doubleK3=0;
doubleK4=0;
y[0]=1;//初值为1
K1=f（x[i],y[i]）;
K2=f（（x[i]+h/2）,（y[i]+h*K1/2））;
K3=f（（x[i]+h/2）,（y[i]+h*K2/2））;
K4=f（（x[i]+h）,（y[i]+h*K3））;
y[i+1]=y[i]+（h/6）*（K1+2*K2+2*K3+K4）;
y1[i]=y[i];
}
for（i=4;i{
y1[i+1]=y[i]+（h/24）*（55*f（x[i],y[i]）-59*f（x[i-1],y[i-1]）+37*f（x[i-2],y[i-2]）-9*f（x[i-3],y[i-3]））;
y[i+1]=y[i]+（h/24）*（9*f（x[i+1],y1[i+1]）+19*f（x[i],y[i]）-5*f（x[i-1],y[i-1]）+f（x[i-2],y[i-2]））;
}
//输出x和y
printf（"\n***********************Adams预报—校正格式求解方程***************************\n"）;
printf（"x[n]y[n]y（精确解）\n\n"）;
for（i=0;i{
printf（"%.4f%.4f%.4f\n",x[i],y[i],y2[i]）;
}
printf（"\n**********************Adams预报—校正格式求解方程****************************\n\n"）;
return0;
}
五、精确解程序如下：
intac_solve（inta,intb,doubleh）
{
doubley[LEN]={0};
doublex[LEN]={0};
inti=0;
doublenum=0;
num=（（b-a）/h）+1;//节点个数
//循环计算出各个节点的坐标
for（i=1;i{
x[0]=a;
x[i]=a+i*h;
}
//带公式计算各个节点的y的值
for（i=0;i{
y[i]=sqrt（x[i]/2+（x[i]*x[i]）/2）;
}
y[0]=1;
//输出x和y
printf（"\n****************************精确解***********************************\n"）;
printf（"\nx[n]y[n]\n\n"）;
for（i=0;i{
printf（"%.4f%.4f\n",x[i],y[i]）;
}
printf（"\n****************************精确解***********************************\n\n"）;
return0;
}
调用前四个函数
#include
#include
#include"Euler.h"
#defineLEN100
intmain（）
{
inta=1;//区间左端点
intb=2;//区间右端点
intcase_flag=0;
doubleh=0.1;//步长
item（）;
scanf（"%d",&case_flag）;
printf（"\n"）;
switch（case_flag）
{
case0:
{
break;
}
case1:
{
Euler（a,b,h）;//Euler法
break;
}
case2:
{
Mid（a,b,h）;//中点格式
break;
}
case3:
{
PreCorr（a,b,h）;//预报—校正格式
break;
}
case4:
{
RungeKutta（a,b,h）;//经典四级四阶R-K格式
break;
}
case5:
{
ac_solve（a,b,h）;//精确解
break;
}
default:
{
printf（"没有您所选择的序号!
!
!
\n"）;
break;
}
}
return0;
}
调用Adams预报—校正格式主函数
intmain（void）
{
inta=1;//区间左端点
intb=2;//区间右端点
doubleh=0.1;//步长
Adams（a,b,h）;//Adams法
return0;
}
五、数值结果及分析
运算结果可用下面表格表示：
表5.1四种方法以及方程精确解的计算结果表
x[i]
Euler格式
中点格式
预报-校正格式
经典四级四阶R-K格式
Adams格式
精确解
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.1000
1.0000
1.1025
1.1025
1.1025
1.1025
1.0747
1.2000
1.0595
1.2100
1.2100
1.2100
1.2100
1.1489
1.3000
1.1210
1.3223
1.3223
1.3223
1.3223
1.2227
1.4000
1.1849
1.4395
1.4395
1.4394
1.4394
1.2961
1.5000
1.2502
1.5613
1.5613
1.5612
1.5612
1.3693
1.6000
1.3166
1.6877
1.6876
1.6876
1.6876
1.4422
1.7000
1.3836
1.8184
1.8184
1.8184
1.8184
1.5149
1.8000
1.4513
1.9535
1.9535
1.9535
1.9535
1.5875
1.9000
1.5194
2.0928
2.0928
2.0927
2.0927
1.6598
2.0000
1.5878
2.2362
2.2361
2.2361
2.2361
1.7321
表5.2四种方法和方程精确解的差值表
x[i]
Euler格式
中点格式
预报-校正格式
经典四级四阶R-K格式
Adams格式
精确解
1.0000
0.0000
0.0000
0.0000
0.0000
0.0000
1.0000
1.1000
0.0747
-0.0278
-0.0278
-0.0278
-0.0278
1.0747
1.2000
0.0894
-0.0611
-0.0611
-0.0611
-0.0611
1.1489
1.3000
0.1017
-0.0996
-0.0996
-0.0996
-0.0996
1.2227
1.4000
0.1112
-0.1434
-0.1434
-0.1433
-0.1433
1.2961
1.5000
0.1191
-0.192
-0.192
-0.1919
-0.1919
1.3693
1.6000
0.1256
-0.2455
-0.2454
-0.2454
-0.2454
1.4422
1.7000
0.1313
-0.3035
-0.3035
-0.3035
-0.3035
1.5149
1.8000
0.1362
-0.366
-0.366
-0.366
-0.366
1.5875
1.9000
0.1404
-0.433
-0.433
-0.4329
-0.4329
1.6598
2.0000
0.1443
-0.5041
-0.504
-0.504
-0.504
1.7321
求和
1.1739
-2.3760
-2.3758
-2.3755
-2.3755
由以上表格可知：
从表5.1和表5.2可以看出在四种方法中Euler格式计算的方程的解与精确解相差最大，即误差较大；中点格式，预报-校正格式和经典四级四阶R-K格式的误差较小，其中经典四级四阶R-K格式离精确解最近，所以，经典四级四阶R-K格式方法在前四种方法中求解方程的解的误差最小，效果最好。
通过计算结果，可知预报格式已是很好的近似，同时Adams内插法又有较高的精度，故我们用相应的外插公式作为预报格式，然后用内插格式进行迭代，从而形成Adams预报—修正格式，所以Adams预报—修正格式效果显著，能达到较高的精度，尽管只迭代一次却与真解十分相近。
第一章习题9实验报告
一、实验目的
通过用Euler法求解高阶常微分方程初值问题
的数值解，深刻领会微分方程数值解的实质，熟练掌握该算法的计算机实现过程，并能从理论及实验结果分析该算法的优缺点。
通过本次实验，学会将一阶常微分方程初值问题的数值解法推广到一阶常微分方程组初值问题中。
这里以Euler法为例进行推广，体会其相同之处及不同之处，同时学会向量值函数的表达。
二、实验内容
用Euler法求解高阶常微分方程初值问题
的数值解，其精确解为
。
三、算法
Euler法求解高阶常微分方程初值问题
其中
为步长。
四、程序
4.1头文件
功能：
头文件中包含了四个方法以及方程精确解求解的具体实现过程，并将各个方法打包成函数，方便在主文件中进行调用。
各个函数分别为
doublef1（doublex,doubley1,doubley2）;
doublef2（doublex,doubley1,doubley2）;
头文件程序：
doublef1（doublex,doubley1,doubley2）
{
returny2;
}
doublef2（doublex,doubley1,doubley2）
{
return-y1;
}
4.2主文件
主文件程序：
intmain（void）
{
doublex[ARRARY_LEN]={0};
doubley[ARRARY_LEN][ARRARY_LEN]={0};
doubley1[ARRARY_LEN][ARRARY_LEN]={0};
inti=0;
inta=0;
intb=1;
doubleh=0.1;
doublenum=0;
num=（（b-a）/h）+1;//节点个数
//循环计算出各个节点的坐标
for（i=1;i{
x[0]=a;
x[i]=a+i*h;
}
//计算精确解
for（i=0;i{
y1[0][0]=1;
y1[1][0]=1;
y1[0][i+1]=sin（x[i+1]）;
y1[1][i+1]=cos（x[i+1]）;
}
//带公式计算各个节点的y的值
for（i=0;i{
y[0][0]=0;
y[1][0]=1;
y[0][i+1]=y[0][i]+h*f1（x[i],y[0][i],y[1][i]）;
y[1][i+1]=y[1][i]+h*f2（x[i],y[0][i],y[1][i]）;
}
//输出x和y
printf（"\n*************************Euler法求解高阶方程*****************************\n"）;
printf（"\nx[n]y1[n]y1（精确解）y2[n]y2（精确解）\n\n"）;
for（i=0;i{
printf（"%6.4f%6.4f%8.4f%6.4f%8.4f\n",x[i],y[0][i],y[1][i],y1[0][i],y1[1][i]）;
}
printf（"\n*************************Euler法求解高阶方程****************************\n\n"）;
return0;
}
五、数值结果及分析
5.1结果截图：
程序运行根据函数调用进行计算:
表5.1Euler法求解高解方程以及方程精确解的计算结果表
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
升级会员 