Matlab解热传导方程代码

SampleMATLABcodes1.%NewtonCoolingLawclear;closeall;clc;h=1;T(1)=10;%T(0)error=1;TOL=1e-6;k=0;dt=1/10;whileerrorTOL,k=k+1;T(k+1)=h*(1-T(k))*dt+T(k);error=abs(T(k+1)-T(k));endt=linspace(0,dt*(k+1),k+1);plot(t,T),holdon,plot(t,1,'r-.')xlabel('Time'),ylabel('Temperature'),title(['T_0=',num2str(T(1)),',T_\infty=1']),legend('CoolingTrend','SteadyState')2.%BoltzmanCoolingLawclear;closeall;clc;h=1;T(1)=10;%T(0)error=1;TOL=1e-6;k=0;dt=1/10000;whileerrorTOL,k=k+1;T(k+1)=h*(1-(T(k))^4)*dt+T(k);error=abs(T(k+1)-T(k));endt=linspace(0,dt*(k+1),k+1);plot(t,T),holdon,plot(t,1,'r-.')xlabel('Time'),ylabel('Temperature'),title(['T_0=',num2str(T(1)),',T_\infty=1']),legend('CoolingTrend','SteadyState')3.%FourierHeatconductionclear;closeall;clc;h=1;n=11;T=ones(n,1);Told=T;T(1)=1;%LeftboundaryT(n)=10;%Rightboundaryx=linspace(0,1,n);dx=x(2)-x(1);dt=dx^2/3;%cflconditionerror=1;TOL=1e-6;k=0;whileerrorTOL,Told=T;k=k+1;fori=2:n-1T(i)=dt*(Told(i+1)-2*Told(i)+Told(i-1))/dx^2+Told(i);enderror=max(abs(T-Told));ifmod(k,5)==0,out(k,:)=T;endendplot(x,out)xlabel('x'),ylabel('Temperature'),title(['FourierHeatConduction']),%legend('CoolingTrend','SteadyState')4.2DHeatEquation%2DHeatEquation.clear;closeall;clcn=10;%gridhasn-2interiorpointsperdimension(overlapping)x=linspace(0,1,n);dx=x(2)-x(1);y=x;dy=dx;TOL=1e-6;T=zeros(n);T(1,1:n)=10;%TOPT(n,1:n)=1;%BOTTOMT(1:n,1)=1;%LEFTT(1:n,n)=1;%RIGHTdt=dx^2/4;error=1;k=0;whileerrorTOLk=k+1;Told=T;fori=2:n-1forj=2:n-1T(i,j)=dt*((Told(i+1,j)-2*Told(i,j)+Told(i-1,j))/dx^2...+(Told(i,j+1)-2*Told(i,j)+Told(i,j-1))/dy^2)...+Told(i,j);endenderror=max(max(abs(Told-T)));endsubplot(2,1,1),contour(x,y,T),title('Temperature(SteadyState)'),xlabel('x'),ylabel('y'),colorbarsubplot(2,1,2),pcolor(x,y,T),shadinginterp,title('Temperature(SteadyState)'),xlabel('x'),ylabel('y'),colorbar5.WaveTranslation%Oscillations-translationleftandrightclear;closeall;clc;forc=[1-1]cc=0;n=261;x=linspace(0,13,n);u=zeros(n,1);u(121:141)=sin(pi*x(121:141));dx=x(2)-x(1);dt=dx;error=1;TOL=1e-6;k=0;whilek110uold=u;k=k+1;fori=2:n-1ifc==1,u(i)=dt*c*(uold(i+1)-uold(i))/dx+uold(i);end%c=1ifc==-1,u(i)=dt*c*(uold(i)-uold(i-1))/dx+uold(i);end%c=-1enderror=max(abs(u-uold));ifmod(k,10)==0,cc=cc+1;out(cc,:)=u;endendifc==1subplot(2,1,1),forhh=1:ccplot(x,out(hh,:)+hh),holdon,endu=zeros(n,1);u(121:141)=sin(pi*x(121:141));plot(x,u)xlabel('u(x)'),ylabel('Time'),title('TranslationtotheLeft')elseifc==-1subplot(2,1,2),forhh=1:ccplot(x,out(hh,:)+hh),holdon,endu=zeros(n,1);u(121:141)=sin(pi*x(121:141));plot(x,u)xlabel('u(x)'),ylabel('Time'),title('TranslationtotheRight')endend6.%waveequationclear;closeall;clc;c=1;n=21;x=linspace(0,1,n);dx=1/(n-1);dt=dx;u(:,1)=sin(pi*x);u(1,2)=0;fori=2:n-1u(i,2)=0.5*(dt^2*c^2*(u(i+1,1)-2*u(i,1)+u(i-1,1))/dx^2+2*u(i,1));endu(n,2)=0;error=1;k=1;whilek100k=k+1;u(1,k+1)=0;fori=2:n-1u(i,k+1)=dt^2*c^2*(u(i+1,k)-2*u(i,k)+u(i-1,k))/dx^2+2*u(i,k)-u(i,k-1);endu(n,k+1)=0;endplot(x,u),xlabel('x'),ylabel('y')