微分方程数值解上机报告常微分方程初值问题数值解法.docx

1、微分方程数值解上机报告常微分方程初值问题数值解法第一章 常微分方程初值问题数值解法实验一1、 实验题目试用(a)欧拉格式(b)中点格式(c)预报校正格式(d)经典四级四阶R-K格式编程计算方程：2、 程序#include#include#includeconst int N=11;double fund(double x,double y);void Euler(double a,double h,double y);void Center(double a,double h,double y);void YuJiao(double a,double h,double y);void SjSj

2、(double a,double h,double y);void Adams(double a,double h,double y);void main()double a,b,h,yN;int i;char option;coutab;couty0;h=(b-a)/(N-1);couta:欧拉法 nb:中心法 nc:预报-校正格式 n;coutd:经典四级四阶R-K格式 n;coute:Adams预报-修正格式 nf:退出n;label: coutoption;switch(option)case a: Euler(a,h,y);break;case b: Center(a,h,y);br

3、eak;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;goto label;cout计算结果为：n xnttynendl;for(i=0;iN;i+)cout a+i*httyiendl;void Euler(double a,double h,double y)int i;for(i=1;iN;i+)yi=yi-1+h*fund(a+(i-1)*h,yi-1);void Cent

4、er(double a,double h,double y)int i;double w;for(i=1;iN;i+)w=fund(a+(i-1)*h,yi-1);yi=yi-1+h*fund(a+(i-1)*h/2,w*h/2+yi-1); void YuJiao(double a,double h,double y)int i;double sun;double y_BeginN=y0;for(i=1;iN;i+)y_Begini=y_Begini-1+h*fund(a+(i-1)*h,y_Begini-1);for(i=1;i0.0001)sun=yi-1+h/2.0*(fund(a+(

5、i-1)*h,yi-1)+fund(a+i*h,y_Begini);y_Begini=sun;yi=sun;void SjSj(double a,double h,double y)int i;double k1,k2,k3,k4;for(i=1;iN;i+)k1=fund(a+(i-1)*h,yi-1);k2=fund(a+(i-1)*h+h/2,yi-1+h*k1/2);k3=fund(a+(i-1)*h+h/2,yi-1+h*k2/2);k4=fund(a+(i-1)*h+h,yi-1+h*k3);yi=yi-1+h/6*(k1+2*k2+2*k3+k4);void Adams(doub

6、le a,double h,double y)int i;double me,xN;double k1,k2,k3,k4;for(i=0;iN;i+)xi=a+i*h;yi=y0;for(i=1;i4;i+)k1=fund(a+(i-1)*h,yi-1);k2=fund(a+(i-1)*h+h/2,yi-1+h*k1/2);k3=fund(a+(i-1)*h+h/2,yi-1+h*k2/2);k4=fund(a+(i-1)*h+h,yi-1+h*k3);yi=yi-1+h/6*(k1+2*k2+2*k3+k4);for(i=4;i0.0001) yi=yi-1+h/24*(9*fund(xi,

7、me)+19*fund(xi-1,yi-1)-5*fund(xi-2,yi-2)+fund(xi-3,yi-3);me=yi;double fund(double x,double y)double s;s=y/(2*x)+x*x/2/y;return s;3、 试验结果及分析(i)运算结果欧拉法解中点法解预报-校正法解经典四级四阶R-K法解1.011111.11.11.100121.10251.10251.21.2051.202831.209971.209971.31.314961.30811.322351.322351.41.42981.41591.439541.439541.51.549

8、41.526181.561411.561411.61.673661.63891.687851.687851.71.802441.754021.818731.818731.81.935621.871511.953931.953931.92.073081.991322.093342.093342.02.21472.113412.236832.23683(ii)结果分析在使用欧拉法数值求解过程中，其计算过程非常简单，即由可直接计算出，由可直接计算出无需用迭代方法求解任何方程，就可求出近似解。通过上面两表的比较，可知随着步长的减小，近似解越来越接逼近精确解，即误差越来越小，当步长充分小时，所得的近似解

9、能足够地逼近精确解。通过中点格式计算出的近似解，格式提高了精度，不过相对于预报-校正格式和经典四级四阶Runge-Kutta格式的高精度，其精度显然还不够高。在上述两表的比较中，可以看出当步长取得适当小，用预报校正格式只需迭代一次就可算出比较好的近似值，得到满足精度要求的解，且比欧拉格式具有更高的精度。经典四级四阶Runge-Kutta格式是高精度的，并且是显示格式，所以该格式具有算法简单并能得到高精度解的优点。通过上表可见，四级四阶Runge-Kutta法的确可计算高精度的解。实验二1、实验题目试用Adams预报修正格式，计算下列方程：3、程序说明：本次实验的程序已嵌入到实验一的程序中、下面

10、是本次实验的主要代码。void Adams(double a,double h,double y)int i;double me,xN;double k1,k2,k3,k4;for(i=0;iN;i+)xi=a+i*h;yi=y0;for(i=1;i4;i+)k1=fund(a+(i-1)*h,yi-1);k2=fund(a+(i-1)*h+h/2,yi-1+h*k1/2);k3=fund(a+(i-1)*h+h/2,yi-1+h*k2/2);k4=fund(a+(i-1)*h+h,yi-1+h*k3);yi=yi-1+h/6*(k1+2*k2+2*k3+k4);for(i=4;i0.0001) yi=yi-1+h/24*(9*fund(xi,me)+19*fund(xi-1,yi-1)-5*fund(xi-2,yi-2)+fund(xi-3,yi-3);me=yi;3、试验结果及分析(i)运行结果11.11.21.31.41.51.61.71.81.92.011.10251.209961.322311.439441.561251.68761.818381.953462.092732.23607(ii)结果分析