使用差分方法求解下面的热传导方程
2(,)4(,)0100.2t xx T x t T x t x t =<<<<
初值条件：2(,0)44T x x x =-
边值条件：(0,)0(1,)0
T t T t ==
使用差分公式
1,,1,222(,)2(,)(,)
2(,)()i j i j i j i j i j i j xx i j T x h t T x t T x h t T T T T x t O h h h -+--++-+=+≈
,1,(,)(,)
(,)()i j i j i j i j
t i j T x t k T x t T T T x t O k k k ++--=+≈
上面两式带入原热传导方程
,1,1,,1,22i j i j
i j i j i j
T T T T T k h +-+--+= 令224k r h
=，化简上式的 ,1,1,1,(12)()i j i j i j i j T r T r T T +-+=-++
i
x j
t j
编程MATLAB 程序，运行结果如下
x t T
function mypdesolution
c=1;
xspan=[0 1];
tspan=[0 0.2];
ngrid=[100 10];
f=@(x)4*x-4*x.^2;
g1=@(t)0;
g2=@(t)0;
[T,x,t]=rechuandao(c,f,g1,g2,xspan,tspan,ngrid);
[x,t]=meshgrid(x,t);
mesh(x,t,T);
xlabel('x')
ylabel('t')
zlabel('T')
function [U,x,t]=rechuandao(c,f,g1,g2,xspan,tspan,ngrid)
% 热传导方程：
% Ut(x,t)=c^2*Uxx(x,t) a<x<b ts<t<tf
% 初值条件：
% u(x,0)=f(x)
% 边值条件：
% u(a,t)=g1(t)
% u(b,t)=g2(t)
%
% 参数说明
% c：方程中的系数
% f：初值条件
% g1,g2：边值条件
% xspan=[a,b]：x的取值范围
% tspan=[ts,tf]：t的取值范围
% ngrid=[n,m]：网格数量，m为x网格点数量，n为t的网格点数量% U：方程的数值解
% x,t：x和t的网格点
n=ngrid(1);
m=ngrid(2);
h=range(xspan)/(m-1);
x=linspace(xspan(1),xspan(2),m);
k=range(tspan)/(n-1);
t=linspace(tspan(1),tspan(2),n);
r=c^2*k/h^2;
if r>0.5
error('为了保证算法的收敛，请增大步长h或减小步长k!')
end
s=1-2*r;
U=zeros(ngrid);
% 边界条件
U(:,1)=g1(t);
U(:,m)=g2(t);
% 初值条件
U(1,:)=f(x);
% 差分计算
for j=2:n
for i=2:m-1
U(j,i)=s*U(j-1,i)+r*(U(j-1,i-1)+U(j-1,i+1));
end
end。