微分方程数值解上机报告常微分方程初值问题数值解法.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(d

2、ouble 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: Cente

3、r(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;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+(

4、i-1)*h,yi-1);void Center(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

5、;i0.0001) sun=yi-1+h/2.0*(fund(a+(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);

6、yi=yi-1+h/6*(k1+2*k2+2*k3+k4); 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*(

7、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; 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

8、.314961.30811.322351.322351.41.42981.41591.439541.439541.51.54941.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法的确可计算高

10、精度的解。实验二1、实验题目试用Adams预报修正格式，计算下列方程：3、程序说明：本次实验的程序已嵌入到实验一的程序中、下面是本次实验的主要代码。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

11、=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)结果分析通过实验一的计算结

12、果，可知预报格式已是很好的近似，同时Adams内插法又有较高的精度，故我们用相应的外插公式作为预报格式，然后用内插格式进行迭代，从而形成Adams预报修正格式，所以Adams预报修正格式效果显著，能达到较高的精度，尽管只迭代一次却与真解十分相近。实验三1、实验题目分别用四级四阶R-K格式法和Adams预报较正格式解高阶方程：2、程序#include#include#includeconst m=2;const n=11;double F1(double x,double y1,double y2);double F2(double x,double y1,double y2);void mai

13、n() double Ymn,Km4; double h,a,b; int i; Y00=0; /初值1 Y10=1; /初值2 a=0.0; /区间左边界 b=1.0; /区间右边界 h=0.1; /步长 for(i=0;in-1;i+) K00=F1(a+i*h,Y0i,Y1i); K10=F2(a+i*h,Y0i,Y1i); K01=F1(a+i*h+h/2,Y0i+h*K00/2,Y1i+h*K10/2); K11=F2(a+i*h+h/2,Y0i+h*K00/2,Y1i+h*K10/2); K02=F1(a+i*h+h/2,Y0i+h*K01/2,Y1i+h*K11/2); K12=

14、F2(a+i*h+h/2,Y0i+h*K01/2,Y1i+h*K11/2); K03=F1(a+i*h+h,Y0i+h*K02,Y1i+h*K12); K13=F2(a+i*h+h,Y0i+h*K02,Y1i+h*K12); Y0i+1=Y0i+h/6*(K00+2*K01+2*K02+K03); Y1i+1=Y1i+h/6*(K10+2*K11+2*K12+K13); cout用R_K法计算结果为：nn; coutXty1tty2n; for(i=0;in;i+) coutx=a+i*ht; coutsetw(8)Y0itY1i; coutn; coutn; cout解析精确解为：nn; c

15、outXty1tty2n; for(i=0;in;i+) coutx=a+i*ht; coutsetw(8)sin(a+i*h)tcos(a+i*h); coutn; double F1(double x,double y1,double y2)return y2;double F2(double x,double y1,double y2)return -y1;3、试验结果及分析(i)运行结果R-K法精确解001010.10.09983330.9950040.09983340.9950040.20.1986680.9800670.1986690.9800670.30.295520.95533

16、70.295520.9553360.40.3894180.9210610.3894180.9210610.50.4794250.8775830.4794260.8775830.60.5646420.8253360.5646420.8253360.70.6442170.7648430.6442180.7648420.80.7173560.6967070.7173560.6967070.90.7833260.6216110.7833270.6216110.841470.5403030.8414710.540302(ii)结果分析经典四级四阶Runge-Kutta格式是高精度的，并且是显示格式，所以

17、该格式具有算法简单并能得到高精度解的优点。通过四级四阶Runge-Kutta格式解高阶方程，由上表可见，四级四阶Runge-Kutta法的确可计算高精度的解。尽管如此，我们仍不能利用计算机求得它们的精确解，这是因为计算机计算，不论运算器能进行多少位数的运算，总会出现舍入误差以及在计算过程中误差的传递。有上表的误差分析一栏，可以看出随着计算的进行，误差在不断的累积。第二章 抛物型方程的差分法实验四1、实验题目用古典显式和古典隐式格式计算抛物型方程 满足初始条件 和边界条件 的近似解，并与解析解进行比较（）。2、程序#include#includeconst double pi=3.1415926

18、;const int N=11;const int M=11;const double t=0.001;const double e=2.718281828459;double Ut(double x);/初始时刻值double Ux1(double time);/左边值double Ux2(double time);/右边值double FUN(double x,double time);void main() int i,k; double UMN; double h,a,b,T1,Tn,r; coutab; coutT1Tn; h=(b-a)/(N-1); r=t/(h*h); for(

19、k=0;kM;k+) Uk0=Ux1(T1+t*k); UkN-1=Ux2(T1+t*k); for(i=0;iN;i+) U0i=Ut(a+h*i); for(k=1;kM;k+) for(i=1;iN-1;i+) Uki=r*Uk-1i-1+(1-2*r)*Uk-1i+r*Uk-1i+1; cout古典显格式计算结果(t=0.005,h=0.1)：n; for(i=0;iN;i+) coutU5i ; coutn; cout精确计算结果(t=0.005,h=0.1)：n; for(i=0;iN;i+) coutFUN(a+i*h,0.005)=1)/为了减小误差当x等于整数时，令x=0；此

20、时sin(x*pi)值不变 x=0; return (sin(pi*x);double Ux1(double time)/左边值return 0;double Ux2(double time)/右边值return 0;double FUN(double x,double time) double s1,s2; if(int)x)=1) x=0; s1=-pow(pi,2)*time; s2=pow(e,s1)*sin(pi*x); return s2;3、试验结果及分析(i)运行结果节点古典显式格式解真解0000.10.2941860.2941380.20.5595750.5594830.30

21、.7701890.7700630.40.9054110.9052630.50.9520050.951850.60.9054110.9052630.70.7701890.7700630.80. 5595750.5594830.90.2941860.294138100(ii)结果分析通过显示差分格式和隐式差分格式求解抛物型方程，由上表的计算结果，可以看出古典隐式格式比古典显式格式具有更高的精度。实验五1、实验题目用Crank-Nicolson差分格式计算抛物型方程 满足初始条件 和边界条件 在处的解，。2、程序#include#includeconst double pi=3.1415926;co

22、nst int N=11;const int M=11;const double t=0.1;const double h=0.1;const double e=2.71828;double Ut(double x);/初始时刻值double Ux1(double time);/左边值double Ux2(double time);/右边值double FUN(double x,double time);void main() int i,k; double U1111,d9; double a,b,T1,Tn,r; double g9,w9; coutab; coutT1Tn; r=t/(h

23、*h); for(k=0;k11;k+) Uk0=Ux1(T1+t*k); Uk10=Ux2(T1+t*k); for(i=0;i11;i+) U0i=Ut(a+h*i); for(k=1;k11;k+) /计算方程常数项 d0=(1-r)*Uk-10+0.5*r*Uk-11; d8=(1-r)*Uk-19+0.5*r*Uk-18; for(i=1;i8;i+) di=(1-r)*Uk-1i+0.5*r*(Uk-1i+1+Uk-1i-1); g0=d0/(1+r); w0=0.5*r/(1+r); for(i=1;i0;i-) Uki=gi-1+wi-1*Uki+1; cout古典显格式计算结果(t=0.1,h=0.1)：n; for(i=0;i11;i+) coutU1i ; coutn; cout精确计算结果(t=0.1,h=0.1)：n; for(i=0;iN;i+) coutFUN(a+i*h,0.1) ; coutn; cout古典显格式计算结果(t=0.2 h=0.1)：n;