姓名学号中国海洋大学偏微分方程的数值解法第八次作业.docx
《姓名学号中国海洋大学偏微分方程的数值解法第八次作业.docx》由会员分享，可在线阅读，更多相关《姓名学号中国海洋大学偏微分方程的数值解法第八次作业.docx（13页珍藏版）》请在冰豆网上搜索。
姓名学号中国海洋大学偏微分方程的数值解法第八次作业
偏微分方程的数值解法
上机习题八
题目：
分别用显格式、隐格式、C-N格式计算方程的近似解，取不同的网比
其中
，
，
，取
，且对显格式，有
。
精确解为
。
求解：
1记差分解为
，精确解为
，
误差为
；并令
，
，
为
维列向量。
2定义矩阵
，
3已知
，给定
，算出相应的
4向前差分格式：
，其中
5向后差分格式：
，其中
6C-N格式：
，其中
程序：
clear;clc;
J=input（'请输入J='）;%50;
h=1/J;%步长
r=input（'请输入网比r='）;%1/6;
tau=r*h^2;%时间步长
T=2;
N=round（T/tau）;
%----------------预设矩阵-----------
I=eye（J-1）;
S=diag（ones（1,J-2）,1）+diag（ones（1,J-2）,-1）;
U0=cos（pi*[h:
h:
1-h]）;U0=U0';
b=zeros（J-1,1）;
F0=ones（J-1,1）;
%---------------------显格式--------------
ifr<=1/2
C1=（1-2*r）*I+r*S;
U1（1:
J-1,1）=U0;
forn=1:
N
b
（1）=exp（-pi^2*n*tau）+1-cos（n*tau）;
b（J-1）=-exp（-pi^2*n*tau）+1-cos（n*tau）;
F=sin（n*tau）*F0;
U1（:
n+1）=C1*U1（:
n）+tau*F+r*b;
end
U1=U1';
end
%---------------------隐格式--------------
C2=inv（（1+2*r）*I-r*S）;
U2（1:
J-1,1）=U0;
forn=2:
N+1
b
（1）=exp（-pi^2*n*tau）+1-cos（n*tau）;
b（J-1）=-exp（-pi^2*n*tau）+1-cos（n*tau）;
F=sin（n*tau）*F0;
U2（:
n）=C2*U2（:
n-1）+tau*C2*F+r*C2*b;
end
U2=U2';
%---------------------C-N格式--------------
C3=inv（（1+r）*I-r/2*S）*（（1-r）*I+r/2*S）;
U3（1:
J-1,1）=U0;
forn=1:
N
b
（1）=exp（-pi^2*（n+1）*tau）+1-cos（（n+1）*tau）+exp（-pi^2*n*tau）+1-cos（n*tau）;
b（J-1）=-exp（-pi^2*（n+1）*tau）+1-cos（（n+1）*tau）-exp（-pi^2*n*tau）+1-cos（n*tau）;
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;
end
U3=U3';
%------------------精确解计算-----------------
x=h:
h:
1-h;
t=0:
tau:
T;
[X,T]=meshgrid（x,t）;
U=exp（-pi^2.*T）.*cos（pi*X）+1-cos（T）;
%-----------------解作图----------------
figure
（1）;
subplot（2,2,1）;
mesh（X,T,U）;title（'精确解图像'）;xlabel（'x'）;ylabel（'t'）;zlabel（'u'）;
k=2;
ifr<=1/2
subplot（2,2,2）;
mesh（X,T,U1）;title（'显格式图像'）;xlabel（'x'）;ylabel（'t'）;zlabel（'u'）;
k=k+1;
end
subplot（2,2,k）;
mesh（X,T,U2）;title（'隐格式图像'）;xlabel（'x'）;ylabel（'t'）;zlabel（'u'）;
subplot（2,2,k+1）;
mesh（X,T,U3）;title（'C-N格式图像'）;xlabel（'x'）;ylabel（'t'）;zlabel（'u'）;
%----------------解误差作图--------------------
figure
（2）;k=1;
ifr<=1/2
k=k+1;
subplot（k,2,1）;mesh（X,T,U1-U）;
title（'显格式误差图像'）;xlabel（'x'）;ylabel（'t'）;zlabel（'u'）;
end
subplot（k,2,k）;mesh（X,T,U2-U）;
title（隐格式误差图像'）;xlabel（'x'）;ylabel（'t'）;zlabel（'u'）;
subplot（k,2,k+1）;mesh（X,T,U3-U）;
title（'C-N格式误差图像'）;xlabel（'x'）;ylabel（'t'）;zlabel（'u'）;
运行结果
（1）网比为1/6
请输入J=40
请输入网比r=1/6
a）解的图像
b）误差图像
（2）网比为1/2
请输入J=40
请输入网比r=1/2
a）解的图像
b）误差图像
（3）网比为1
请输入J=40
请输入网比r=1
a）解的图像
b）误差图像
（4）网比为2
请输入J=40
请输入网比r=2
a）解的图像
b）误差图像
（5）网比为5
请输入J=40
请输入网比r=5
a）解的图像
b）误差图像
结果分析
1从解的图像和误差图像可以看出，在
附近，解的误差较大。
2网比从
变化到
，误差增大。
数量级从
增大到
。
升级会员 