用Crank-Nicolson差分格式计算抛物型方程.docVIP

实验二（习题2.2） 题目 用Crank-Nicolson差分格式计算抛物型方程 满足初始条件 和边界条件 在处的解，。 程序 #includeiostream.h #includemath.h const 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 U[11][11],d[9]; double a,b,T1,Tn,r; double g[9],w[9]; cout请输入x所属区间[a,b]\n; cinab; cout请输入t所属区间(t1,tn)\n; cinT1Tn; r=t/(h*h); for(k=0;k11;k++) { U[k][0]=Ux1(T1+t*k); U[k][10]=Ux2(T1+t*k); } for(i=0;i11;i++) U[0][i]=Ut(a+h*i); for(k=1;k11;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;i8;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;i9;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;i0;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;i11;i++) { coutU[1][i] ; } cout\n; cout精确计算结果(t=0.1,h=0.1)：\n; for(i=0;iN;i++) { coutFUN(a+i*h,0.1) ; } cout\n; cout古典显格式计算结果(t=0.2 h=0.1)：\n; for(i=0;i11;i++) { coutU[2][i] ; } cout\n; cout精确计算结果(t=0.2,h=0.1)：\n; for(i=0;iN;i++) { coutFUN(a+i*h,0.2) ; } cout\n; } 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; } 运行结果 节点 t=0.1 t=0.2 数值解 真解 数值解 真解 0 0 0 0 0 0.1 0.264775 0.115173 0.12841 0.042926 0.2 0.273488 0.219072 0.0177272 0.08165 0.3 0.305345 0.301527 0.113697 0.112382 0.4 0.338249 0.354466 0.154565 0.132113 0.5 0.356193 0.372708 0.164228 0.13891