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

第一章 常微分方程初值问题数值解法实验一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(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);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+(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;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); 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*(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.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)结果分析在使用欧拉法数值求解过程中，其计算过程非常简单，即由可直接计算出，由可直接计算出无需用迭代方法求解任何方程，就可求出近似解。通过上面两表的比较，可知随着步长的减小，近似解越来越接逼近精确解，即误差越来越小，当步长充分小时，所得的近似解能足够地逼近精确解。通过中点格式计算出的近似解，格式提高了精度，不过相对于预报-校正格式和经典四级四阶Runge-Kutta格式的高精度，其精度显然还不够高。在上述两表的比较中，可以看出当步长取得适当小，用预报校正格式只需迭代一次就可算出比较好的近似值，得到满足精度要求的解，且比欧拉格式具有更高的精度。经典四级四阶Runge-Kutta格式是高精度的，并且是显示格式，所以该格式具有算法简单并能得到高精度解的优点。通过上表可见，四级四阶Runge-Kutta法的确可计算高精度的解。实验二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=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)结果分析通过实验一的计算结果，可知预报格式已是很好的近似，同时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 main()double Ymn,Km4;double h,a,b;int i;Y00=0; /初值1Y10=1; /初值2a=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=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; coutXty1tty2n;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.9553370.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格式是高精度的，并且是显示格式，所以该格式具有算法简单并能得到高精度解的优点。通过四级四阶Runge-Kutta格式解高阶方程，由上表可见，四级四阶Runge-Kutta法的确可计算高精度的解。尽管如此，我们仍不能利用计算机求得它们的精确解，这是因为计算机计算，不论运算器能进行多少位数的运算，总会出现舍入误差以及在计算过程中误差的传递。有上表的误差分析一栏，可以看出随着计算的进行，误差在不断的累积。第二章 抛物型方程的差分法实验四1、实验题目用古典显式和古典隐式格式计算抛物型方程 满足初始条件 和边界条件 的近似解，并与解析解进行比较（）。2、程序#include#includeconst double pi=3.1415926;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(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；此时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.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;const 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*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;for(i=0;i11;i+) coutU2i ; coutn;cout精确计算结果(t=0.2,h=0.1)：n; for(i=0;iN;i+) coutFUN(a+i*h,0.2) ; coutn;double Ut(double x)/初始时刻值if(x-int(x)=0)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(x-int(x)=0)x=0; s1=-pow(pi,2)*time; s2=pow(e,s1)*sin(pi*x); return s2;3、试验结果及分析(i)运行结果节点t=0.1t=0.2数值解真解数值解真解000000.10.2647750.1151730.128410.0429260.20.2734880.2190720.01772720.081650.30.3053450.3015270.1136970.1123820.40.3382490.3544660.1545650.1321130.50.3561930.3727080.1642280.1389110.60.3482590.3544660.1540490.1321130.70.3078640.3015270.1293180.1123820.80.2319270.2190720.09325980.081650.90.1197650.1151730.04982320.04292610000（古典显格式计算结果(t=0.2,h=0.1)和精确计算结果(t=0.2,h=0.1)(ii)结果分析通过上表的计算结果，可以看出用Crank-Nicolson差分格式(三对角追赶法)计算得出的结果误差比较大，精度并不十分理想，但随着步长的缩小，结果精度会有所提高。第三章 椭圆形方程的差分方法实验六1、实验题目利用五点差分格式近似Dirichlet问题 取步长1，试用Jacobi迭代、Guass-Seidel迭代和SOR迭代求解。2、程序#includecmath#includeiostream#includeiomanipconst int N=7;const int M=9;const double pi=3.1415926;double UpFy(double x,double y);double DoFy(double x,double y);double LeFx(double x,double y);double RiFx(double x,double y);using namespace std;int main() static double Matrix_1NM,Matrix_2NM;double x1,x2,y1,y2;/区域边界double dx,dy;/步长double Max,e=1.0;/判断标准(误差)int i,j;cout请输入矩形区域区间:n;cout: ;cinx1x2;cout: ;ciny1y2; dx=(x2-x1)/(M-1);dy=(y2-y1)/(N-1);for(i=0;iN;i+) Matrix_1i0=LeFx(x1,i*dy+y1);Matrix_1iM-1=RiFx(x2,i*dy+y1); for(j=0;j0.00001) for(i=1;iN-1;i+)for(j=1;jM-1;j+)Matrix_2ij=(Matrix_1i-1j+Matrix_1i+1j+Matrix_1ij+1+Matrix_1ij-1)/4.0; Max=fabs(Matrix_111-Matrix_211);/求出误差项最大的元素作为判断标准 for(i=1;iN-1;i+)for(j=1;jMax) Max=fabs(Matrix_2ij-Matrix_1ij); e=Max; for(i=1;iN-1;i+)/为下一次循环做准备for(j=1;jM-1;j+) Matrix_1ij=Matrix_2ij; for(i=0;iN;i+) for(j=0;jM;j+)coutsetw(8)setiosflags(ios:left)Matrix_1ij ; cout1)/当x为4的整数倍时为了提高精度，令x=0x=0;return (sin(x/4.0*pi);double LeFx(double x,double y)return (y*(3-y);double RiFx(double x,double y)return 0;3、试验结果及分析(h=0.5,k=0.5)(i)运行结果Y=00 0.382683 0.707107 0.92388 1 0.92388 0.707107 0.382683 0Y=0.51.25 0.975376 0.896633 0.869842 0.81858 0.70552 0.521264 0.27740Y=1.02 1.37219 1.03422 0.840292 0.698975 0.558373 0.395042 0.20566 0Y=1.52.25 1.4792 1.02777 0.758162 0.578685 0.433982 0.294893 0.150210Y=2.02 1.26683 0.839537 0.58593 0.423655 0.304011 0.200362 0.10030Y=2.51.25 0.748615 0.477633 0.322393 0.226026 0.158071 0.102265 0.0506396 0Y=3.00 0 0 0 0 0 0 0 0(ii)结果分析由运行结果可以看出用Jacobi迭代法计算的结果可以达到预期的精度，从而说明 Jacobi迭代在运算速度和精度上都比较好，是解决这一类问题一个不错的数值方法。第四章 双曲型方程的差分方法实验七1、 实验题目给出初边值问题 用显式差分格式求数值解。2、 程序#include#includeusing namespace std;const int N=11;const int M=11;const double pi=3.1415926;double StartValue(double x);double RightValue(double t);double LeftValue(double t);void main()static double MQMN,MZMN,MJMN;double h,t;int i,j;h=0.1;t=0.1;double R=(t/h)*(t/h);for(i=0;iM;i+) MQi0=LeftValue(i*t); MQiN-1=RightValue(i*t); MQi0=LeftValue(i*t); MQiN-1=RightValue(i*t); MQi0=LeftValue(i*t); MQiN-1=RightValue(i*t); for(i=0;iN;i+)MZ0i=StartValue(i*h); for(i=1;iN-1;i+) MZ1i=R/2.0*(MZ0i-1+MZ0i+1)+(1-R)*MZ0i; for(i=2;iM;i+)for(j=1;jN-1;j+)MZij=R*MZi-1j+1+2*(1-R)*MZi-1j+MZi-1j-1-MZi-2j;for(i=0;iN;i+)MQ0i=StartValue(i*h); MQ1i=MQ0i;for(i=2;iM;i+)for(j=1;jN-1;j+)MQij=R*MQi-1j+1+2*(1-R)*MQi-1j+MQi-1j-1-MQi-2j; for(i=0;iM;i+)for(j=0;jN;j+)coutMQij ;coutnendl; for(i=0;iM;i+)for(j=0;jN;j+)MJij=1/8.0*cos(pi*i*t)*sin(pi*j*h); cout中心结果：n;for(i=0;iM;i+)for(j=0;jN;j+)coutMZij ;coutnendl; cout精确解n; for(i=0;iM;i+)for(j=0;jN;j+)coutMJij ;coutnendl;/初值问题函数double StartValue(double x)if(int)(x)/x=1)x=0; double result=1/8.0*sin(pi*x);return result;/左边值问题函数double LeftValue(double t)return 0;/右边值问题函数double RightValue(double t)return 0;3、 试验结果及分析(i)运算结果u 初边值使用向前差商法T=00 0.0386271 0.0734732 0.101127 0.118882 0.125 0.118882 0.101127 0.0734732 0.0386271 0T=0.10 0.0386271 0.0734732 0.101127 0.118882 0.125 0.118882 0.101127 0.0734732 0.0386271 0T=0.20 0.034846 0.0662811 0.0912281 0.107245 0.112764 0.107245 0.0912281 0.0662811 0.034846 0T=0.30 0.027654 0.052601 0.072399 0.0851102 0.0894901 0.0851102 0.072399 0.052601 0.027654 0T=0.40 0.0177549 0.0337719 0.046483 0.0546441 0.0574562 0.0546441 0.046483 0.0337719 0.0177549 0T=0.50 0.00611794 0.011637 0.016017 0.0188291 0.0197981 0.0188291 0.016017 0.011637 0.00611793 0T=0.60 -0.00611793 -0.011637 -0.016017 -0.0188291 -0.0197981 -0.0188291 -0.016017 -0.011637 -0.00611794 0T=0.70 -0.0177549 -0.0337719 -0.046483 -0.0546441 -0.0574562 -0.0546441 -0.046483 -0.0337719 -0.0177549 0T=0.80 -0.027654 -0.052601 -0.072399 -0.0851102 -0.0894901 -0.0851102 -0.072399 -0.052601 -0.027654 0T=0.90 -0.034846 -0.0662811 -0.0912281 -0.107245 -0.112764 -0.107245 -0.0912281 -0.0662811 -0.034846 0T=1.00 -0.0386271 -0.0734732 -0.101127 -0.118882 -0.125 -0.118882 -0.101127 -0.0734732 -0.0386271 0u 初边值使用中心差商法结果：T=0.0 0.0386271 0.0734732 0.101127 0.118882 0.125 0.118882 0.101127 0.0734732 0.0386271 0T=0.10 0.0367366 0.0698771 0.0961776 0.113064 0.118882 0.113064 0.0961776 0.0698771 0.0367366 0T=0.20 0.03125 0.059441 0.0818136 0.0961776 0.101127 0.0961776 0.0818136 0.059441 0.03125 0T=0.30 0.0227045 0.0431864 0.059441 0.0698771 0.0734732 0.0698771 0.059441 0.0431864 0.0227045 0T=0.40 0.0119364 0.0227045 0.03125 0.0367366 0.0386271 0.0367366 0.03125 0.0227045 0.0119364 0T=0.50 1.03501e-009 1.96871e-009 2.70969e-009 3.18543e-009 -1.38778e-017 -3.18543e-009 -2.70969e-009 -1.96871e-009 -1.03501e-009 0T=0.60 -0.0119364 -0.0227045 -0.03125 -0.0367366 -0.0386271 -0.0367366 -0.03125 -0.0227045 -0.0119364 0T=0.70 -0.0227045 -0.0431864 -0.059441 -0.0698771 -0.0734732 -0.0698771 -0.059441 -0.0431864 -0.0227045 0T=0.80 -0.03125 -0.059441 -0.0818136 -0.0961776 -0.101127 -0.0961776 -0.0818136 -0.059441 -0.03125 0T=0.90 -0.0367366 -0.0698771 -0.0961776 -0.113064 -0.118882 -0.113064 -0.0961776 -0.0698771 -0.0367366 0T=1.00