热传导方程向后差分格式的MATLAB程序

向后差分格式MATLAB编程：

c lear;clc;

format short e

a=input('请输入系数a的值');

l=input('请输入长度l的值');

M=input('请输入将区间[0,1]等分的个数M ');

ot=input('请输入时间增量ot的值');

n=input('请输入运行次数n的值');

ox=1/M;

x0=zeros(M+1,1)

for ii=1:M

x0(ii+1)=ii*ox;

end

u=sin(pi*x0/l);

r=a*ot/(ox)^2;

for ii=1:n

%数据的输入

B=zeros(M-1,1);

A=zeros(M-2,1);

C=zeros(M-2,1);

S=zeros(M-1,1);

for ii=1:M-2

B(ii)=1+2*r;A(ii)=-r;C(ii)=-r;

S(ii)=u(ii+1,1);

end

B(M-1,1)=1+2*r;S(M-1,1)=u(M,1);u(1,2)=0;u(M+1,2)=0;

S(1,1)=S(1,1)+r*u(1,2);S(M-1,1)=S(M-1,1)+r*u(M+1,2);

%追赶法

S(1)=S(1)/B(1);T=B(1);k=2;

while k~=M

B(k-1)=C(k-1)/T;

T=B(k)-A(k-1)*B(k-1);

S(k)=(S(k)-A(k-1)*S(k-1))/T;

k=k+1

end

k=1;

while k~=M-1

S(M-1-k)=S(M-1-k)-B(M-1-k)*S(M-k);

k=k+1;

end

u(2:M,2)=S;

u(:,1)=u(:,2);

end

%计算精确解

for x=0:M

u(x+1,2)=exp(-(pi*a/l)^2*n*ot)*sin(pi*x*ox/l); end

fprintf('最后时刻数值解与精确解分别为');disp(u);