巢湖学院2013届本科毕业论文（设计）
附页
算例2.1.1的编程： a?2
a?111取h?，??，则r?2?，满足稳定性条件
h21010a?111另取h?，??，则r?2?，亦满足稳定性条件
h22020a?111另取h?，??，则r?2?，亦满足稳定性条件
h24040format long a=2; l=1; T=1; N=10; M=10; h=l/N; to=T/M;
r=(a*to)/h^2; for j=1:N+1
x(j)=(j-1)*h; for k=1:M+1
t(k)=(k-1)*to;
u(j,k)=exp(x(j)+2*t(k)); end end
u %求解精确解
for j=1:N+1
x(j)=(j-1)*h;
us(j,1)=exp(x(j)); end
for k=1:M+1
15
二维双曲型方程的交替方向隐格式解法
t(k)=(k-1)*to;
us(1,k)=exp(2*t(k));
us(N+1,k)=exp(1+2*t(k)); end
for k=2:M+1 for j=2:N
us(j,k)=r*us(j-1,k-1)+(1-2*r)*u(j,k-1)+r*us(j+1,k-1); end end
us %求解数值解
for k=1:M+1 for j=1:N+1
R(j,k)=abs(u(j,k)-us(j,k)); end end
R %计算误差
Rmax=max(max(R)) %求误差的最大值
精确解与数值解的比较： x=0:0.1:1; hold on
plot(x,u(:,M+1),'b'); plot(x,us(:,M+1),'y');
title('t=1，h=1/10，τ=1/10时精确解和数值解的比较') text(0.05,21,'蓝：精确解'); text(0.05,20,'黄：数值解'); hold off
取不同步长时的误差比较： x=0:1/10:1; y=0:1/20:1; z=0:1/40:1; hold on
plot(x,R(:,M+1),'b'); hold off
M分别取10，20，40
16