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

向后差分格式MATLAB编程：clear;clc;formatshortea=input('请输入系数a的值');l=input('请输入长度l的值');M=input('请输入将区间[0,1]等分的个数M');ot=input('请输入时间增量ot的值');n=input('请输入运行次数n的值');ox=1/M;x0=zeros(M+1,1)forii=1:Mx0(ii+1)=ii*ox;endu=sin(pi*x0/l);r=a*ot/(ox)^2;forii=1:n%数据的输入B=zeros(M-1,1);A=zeros(M-2,1);C=zeros(M-2,1);S=zeros(M-1,1);forii=1:M-2B(ii)=1+2*r;A(ii)=-r;C(ii)=-r;S(ii)=u(ii+1,1);endB(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;whilek~=MB(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+1endk=1;whilek~=M-1S(M-1-k)=S(M-1-k)-B(M-1-k)*S(M-k);k=k+1;endu(2:M,2)=S;u(:,1)=u(:,2);end%计算精确解forx=0:Mu(x+1,2)=exp(-(pi*a/l)^2*n*ot)*sin(pi*x*ox/l);endfprintf('最后时刻数值解与精确解分别为');disp(u);