《姓名学号中国海洋偏微分方程的数值解法第八次作业》由会员分享，可在线阅读，更多相关《姓名学号中国海洋偏微分方程的数值解法第八次作业（11页珍藏版）》请在金锄头文库上搜索。
1、偏微分方程的数值解法偏微分方程的数值解法上机习题八上机习题八题目：题目：分别用显格式、隐格式、C-N格式计算方程的近似解，取不同的网比 r 2222sin ,01,02,0cos0,1 cos1,1 costtUUtxttx U xxUtetUtet 其中 ， ， ，取，且对显格式，有 。2rh1hJntn1 1=1,2,56 2r，1 2r 精确解为 。2,cos1 costU x text 求解：求解：1记差分解为，精确解为，121,Tnnnn JUUUUL121,Tnnnn JuuuuL误差为；并令，nnnEUu121,= sin,sin,sinTTnnnn JnnnFffftttLL，
2、 为 维220,0,0,=1 cos,0,0,1 cosnnTTttnnn Jnnbuuetet LLnb1J 列向量。2定义矩阵 ，1 1 1I 01 101 10S 3已知，给定，算出相应的 2T ,J r, ,hN4向前差分格式： ，其中111,0,1,1nnnnUCUA FrA bnNK1,1 2CA B AI Br IrS5向后差分格式：，其中-111,1,nnnnUCUA FrA bnN K1,1+2,CA B Ar IrS BI6C-N格式：，+1-11111,0,1,122nnnnnnrUCUAFFAbbnNK其中1,1+,122rrCA B Ar IS Br IS程序：程序：
3、clear;clc;clear;clc; J=input(J=input( 请输入请输入J=J=););%50;%50; h=1/J;h=1/J;% % 步长步长 r=input(r=input( 请输入网比请输入网比r=r=););%1/6;%1/6; tau=r*h2;tau=r*h2;% % 时间步长时间步长T=2;T=2; N=round(T/tau);N=round(T/tau); %-%-预设矩阵预设矩阵-I=eye(J-1);I=eye(J-1); S=diag(ones(1,J-2),1)+diag(ones(1,J-2),-1);S=diag(ones(1,J-2),1)+d
4、iag(ones(1,J-2),-1); U0=cos(pi*h:h:1-h);U0=U0;U0=cos(pi*h:h:1-h);U0=U0; b=zeros(J-1,1);b=zeros(J-1,1); F0=ones(J-1,1);F0=ones(J-1,1); %-%-显格式显格式-ifif r=1/2r=1/2C1=(1-2*r)*I+r*S;C1=(1-2*r)*I+r*S;U1(1:J-1,1)=U0;U1(1:J-1,1)=U0;forfor n=1:Nn=1:Nb(1)=exp(-pi2*n*tau)+1-cos(n*tau);b(1)=exp(-pi2*n*tau)+1-co
5、s(n*tau);b(J-1)=-exp(-pi2*n*tau)+1-cos(n*tau);b(J-1)=-exp(-pi2*n*tau)+1-cos(n*tau);F=sin(n*tau)*F0;F=sin(n*tau)*F0;U1(:,n+1)=C1*U1(:,n)+tau*F+r*b;U1(:,n+1)=C1*U1(:,n)+tau*F+r*b;endendU1=U1;U1=U1; endend %-%-隐格式隐格式-C2=inv(1+2*r)*I-r*S);C2=inv(1+2*r)*I-r*S); U2(1:J-1,1)=U0;U2(1:J-1,1)=U0; forfor n=2:N
6、+1n=2:N+1b(1)=exp(-pi2*n*tau)+1-cos(n*tau);b(1)=exp(-pi2*n*tau)+1-cos(n*tau);b(J-1)=-exp(-pi2*n*tau)+1-cos(n*tau);b(J-1)=-exp(-pi2*n*tau)+1-cos(n*tau);F=sin(n*tau)*F0;F=sin(n*tau)*F0;U2(:,n)=C2*U2(:,n-1)+tau*C2*F+r*C2*b;U2(:,n)=C2*U2(:,n-1)+tau*C2*F+r*C2*b; endend U2=U2;U2=U2; %-C-N%-C-N格式格式-C3=inv(
7、1+r)*I-r/2*S)*(1-r)*I+r/2*S);C3=inv(1+r)*I-r/2*S)*(1-r)*I+r/2*S); U3(1:J-1,1)=U0;U3(1:J-1,1)=U0; forfor n=1:Nn=1:Nb(1)=exp(-pi2*(n+1)*tau)+1-cos(n+1)*tau)+exp(-pi2*n*tau)+1-cos(n*tau);b(1)=exp(-pi2*(n+1)*tau)+1-cos(n+1)*tau)+exp(-pi2*n*tau)+1-cos(n*tau);b(J-1)=-exp(-pi2*(n+1)*tau)+1-cos(n+1)*tau)-ex
8、p(-pi2*n*tau)+1-cos(n*tau);b(J-1)=-exp(-pi2*(n+1)*tau)+1-cos(n+1)*tau)-exp(-pi2*n*tau)+1-cos(n*tau);F=(sin(n*tau)+sin(n+1)*tau)*F0;F=(sin(n*tau)+sin(n+1)*tau)*F0;U3(:,n+1)=C3*U3(:,n)+tau/2*inv(1+r)*I-r/2*S)*F+r/2*inv(1+r)*I-r/2*S)*b;U3(:,n+1)=C3*U3(:,n)+tau/2*inv(1+r)*I-r/2*S)*F+r/2*inv(1+r)*I-r/2*S
9、)*b;endend U3=U3;U3=U3; %-%-精确解计算精确解计算-x=h:h:1-h;x=h:h:1-h; t=0:tau:T;t=0:tau:T; X,T=meshgrid(x,t);X,T=meshgrid(x,t); U=exp(-pi2.*T).*cos(pi*X)+1-cos(T);U=exp(-pi2.*T).*cos(pi*X)+1-cos(T); %-%-解作图解作图-figure(1);figure(1); subplot(2,2,1);subplot(2,2,1); mesh(X,T,U);title(mesh(X,T,U);title( 精确解图像精确解图像
10、);xlabel();xlabel(xx);ylabel();ylabel(tt);zlabel();zlabel(uu););k=2;k=2; ifif r=1/2r=1/2 subplot(2,2,2);subplot(2,2,2); mesh(X,T,U1);title(mesh(X,T,U1);title( 显格式图像显格式图像 );xlabel();xlabel(xx);ylabel();ylabel(tt);zlabel();zlabel(uu););k=k+1;k=k+1; endend subplot(2,2,k);subplot(2,2,k); mesh(X,T,U2);ti
11、tle(mesh(X,T,U2);title( 隐格式图像隐格式图像 );xlabel();xlabel(xx);ylabel();ylabel(tt);zlabel();zlabel(uu););subplot(2,2,k+1);subplot(2,2,k+1); mesh(X,T,U3);title(mesh(X,T,U3);title(C-NC-N格式图像格式图像 );xlabel();xlabel(xx);ylabel();ylabel(tt);zlabel();zlabel(uu);); %-%-解误差作图解误差作图-figure(2);k=1;figure(2);k=1; ifif
12、 r=1/2r=1/2k=k+1;k=k+1; subplot(k,2,1);mesh(X,T,U1-U);subplot(k,2,1);mesh(X,T,U1-U); title(title( 显格式误差图像显格式误差图像 );xlabel();xlabel(xx);ylabel();ylabel(tt);zlabel();zlabel(uu););endend subplot(k,2,k);mesh(X,T,U2-U);subplot(k,2,k);mesh(X,T,U2-U);title(title(隐格式误差图像隐格式误差图像 );xlabel();xlabel(xx);ylabel(
13、);ylabel(tt);zlabel();zlabel(uu););subplot(k,2,k+1);mesh(X,T,U3-U);subplot(k,2,k+1);mesh(X,T,U3-U); title(title(C-NC-N格式误差图像格式误差图像 );xlabel();xlabel(xx);ylabel();ylabel(tt);zlabel();zlabel(uu););运行结果运行结果（1） 网比为网比为1/6请输入J=40 请输入网比r=1/6a) 解的图像解的图像b) 误差图像误差图像（2） 网比为网比为1/2请输入J=40 请输入网比r=1/2a) 解的图像解的图像b) 误差图像误差图像（3） 网比为网比为1请输入J=40 请输入网比r=1a) 解的图像解的图像b) 误差图像误差图像（4） 网比为网比为2请输入J=40 请输入网比r=2a) 解的图像解的图像b) 误差图像误差图像（5） 网比为网比为5请输入J=40 请输入网比r=5a) 解的图像解的图像b) 误差图像误差图像结果分析结果分析1从解的图像和误差图像可以看出，在附近，解的误差较大。0t2网比从变化到，误差增大。数量级从增大到。615310210
