偏微分方程向后差分问题（附图）

moller

扩散方程系数为1
古典隐式向后差分格式
代码如下
精确解exp(-X+T);
c1=exp(t);% c1=u(0,t)
c2=exp(-1+t);% c2=u(a,t)
a=1;% 0<x<a
b=1;% 0<t<b
c=1;% 常系数为1
m=5001;% t轴节点数
n=51;% x轴节点数

达人看看哪里出错了，图像在后~~
本人非数学专业，对偏微分的数值解理解不深
希望有人指点迷津，不慎感激！！

% 向后差分格式

clear;
% 初始化
a=1;% 0<x<a
b=1;% 0<t<b
c=1;% 常系数为1
m=5001;% t轴节点数
n=51;% x轴节点数
h=a/(n-1);
k=b/(m-1);% 分别计算步长h，时间步长tao(k)
r=c^2*k/h^2;% 计算网格比
s=1+2*r;% 系数
x=[0:h:1];
t=[0:k:1];
[X,T]=meshgrid(x,t);
c1=exp(t);% c1=u(0,t)
c2=exp(-1+t);% c2=u(a,t)
U=zeros(n,m);
q=zeros(n-1,1);

% 边界条件
U(1,1:m)=c1;
U(n,1:m)=c2;

% 矩阵A
e=ones(n-1,1);
A=spdiags([-r*e,(1+2*r)*e,-r*e],[-1,0,1],n-1,n-1);

U(2:n-1,1)=exp(-(h:h:(n-2)*h))';
for j=2:m
    % 向量q
    q(1)=U(1,j-1)+r*c1(j);
    q(n-1)=U(n-1,j-1)+r*c2(j);
    q(2:n-2)=U(2:n-2,j-1);
    U2=A\q;
    U3(:,j)=U2;   
end
U=[U3;c2];
U=U';

% 画图
%figure(1);
subplot(221);
mesh(X,T,U);
title('图1');
subplot(222);
%figure(2);
U1=exp(-X+T);
mesh(X,T,U1);
title('图2');
%figure(3);
subplot(212);
mesh(X,T,abs(U1-U));
title('图3');
legend('向后差分[n=51]','','图1：数值求解','','图2：精确解','','图3：绝对误差',-1);

[ 本帖最后由 eight 于 2007-9-15 10:26 编辑 ]

moller

从图片来看
是哪个环节出问题了呢？？

moller

顶起来~~

moller

还是没有人理啊~~

xjzuo

若没有记错的话, 这个问题似乎在别的论坛已经解决过了.