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

向前差分格式MATLAB编程： clear;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

D=(1-2*r)*eye(M-1);

temp=r*linspace(1,1,M-2);

D=D+diag(temp,1)+diag(temp,-1); S=D*S

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);