导热方程求解matlab

使用差分方法求解下面的热传导方程2(,)4(,)0100.2txxTxtTxtxt初值条件：2(,0)44Txxx边值条件：(0,)0(1,)0TtTt使用差分公式1,,1,222(,)2(,)(,)2(,)()ijijijijijijxxijTxhtTxtTxhtTTTTxtOhhh,1,(,)(,)(,)()ijijijijtijTxtkTxtTTTxtOkkk上面两式带入原热传导方程,1,1,,1,22ijijijijijTTTTTkh令224krh，化简上式的,1,1,1,(12)()ijijijijTrTrTTixjt1,ijrT,(12)ijrT1,ijrT,1ijT编程MATLAB程序，运行结果如下00.20.40.60.8100.050.10.150.200.20.40.60.81xtTfunctionmypdesolutionc=1;xspan=[01];tspan=[00.2];ngrid=[10010];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)axbtsttf%初值条件：%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;ifr0.5error('为了保证算法的收敛，请增大步长h或减小步长k!')ends=1-2*r;U=zeros(ngrid);%边界条件U(:,1)=g1(t);U(:,m)=g2(t);%初值条件U(1,:)=f(x);%差分计算forj=2:nfori=2:m-1U(j,i)=s*U(j-1,i)+r*(U(j-1,i-1)+U(j-1,i+1));endend