MATLAB编辑一维热传导方程的模拟程序

求解下列热传导问题：

2T1T0zL02tz2 Tz,01z

T0,t1,TL,t01L1,

程序：

function heat_conduction() %一维齐次热传导方程

options={'空间杆长L','空间点数N' ,'时间点数M','扩散系数alfa','稳定条件的值lambda(取值必须小于0.5)',};

topic='seting';

lines=1;

def={'1','100','1000','1','0.5'};

h=inputdlg(options,topic,lines,def);

L=eval(h{1});

N=eval(h{2});

M=eval(h{3});

alfa=eval(h{4});

lambda=eval(h{5});%lambda的值必须小于0.5

%***************************************************

h=L/N;%空间步长

z=0:h:L;

z=z';

tao=lambda*h^2/alfa;%时间步长

tm=M*tao;%热传导的总时间tm

t=0:tao:tm;

t=t';

%计算初值和边值

T=zeros(N+1,M+1);

Ti=init_fun(z);

To=border_funo(t);

Te=border_fune(t);

T(:,1)=Ti;

T(1,:)=To;

T(N+1,:)=Te;

%用差分法求出温度T与杆长L、时间t的关系

for k=1:M

m=2;

while m

T(m,k+1)=lambda*(T(m+1,k)+T(m-1,k))+(-2*lambda+1)*T(m,k); m=m+1;

end;

end;

%设置立体网格

for i=1:M+1

X(:,i)=z;

end;

for j=1:N+1

Y(j,:)=t;

end

mesh(X,Y,T);

view([1 -1 1]);

xlabel('Z');

ylabel('t');

zlabel('T');

function y=init_fun(z)%初值条件

y=1-z.^2;

return

function y=border_funo(t)%z=0的边界条件

y=1+t.*0;

return

function y=border_fune(t)%z=L的边界条件

y=t*.0;

return

运行情况：

按“run”运行时，弹出窗口

将图中相关数据更改为：

点击图框中的“OK”，在“command window”中输出结果为：

求解下列热传导问题：

2T1T0zL02tz2 Tz,01z

T0,t1,TL,t01L1,

程序：

function heat_conduction() %一维齐次热传导方程

options={'空间杆长L','空间点数N' ,'时间点数M','扩散系数alfa','稳定条件的值lambda(取值必须小于0.5)',};

topic='seting';

lines=1;

def={'1','100','1000','1','0.5'};

h=inputdlg(options,topic,lines,def);

L=eval(h{1});

N=eval(h{2});

M=eval(h{3});

alfa=eval(h{4});

lambda=eval(h{5});%lambda的值必须小于0.5

%***************************************************

h=L/N;%空间步长

z=0:h:L;

z=z';

tao=lambda*h^2/alfa;%时间步长

tm=M*tao;%热传导的总时间tm

t=0:tao:tm;

t=t';

%计算初值和边值

T=zeros(N+1,M+1);

Ti=init_fun(z);

To=border_funo(t);

Te=border_fune(t);

T(:,1)=Ti;

T(1,:)=To;

T(N+1,:)=Te;

%用差分法求出温度T与杆长L、时间t的关系

for k=1:M

m=2;

while m

T(m,k+1)=lambda*(T(m+1,k)+T(m-1,k))+(-2*lambda+1)*T(m,k); m=m+1;

end;

end;

%设置立体网格

for i=1:M+1

X(:,i)=z;

end;

for j=1:N+1

Y(j,:)=t;

end

mesh(X,Y,T);

view([1 -1 1]);

xlabel('Z');

ylabel('t');

zlabel('T');

function y=init_fun(z)%初值条件

y=1-z.^2;

return

function y=border_funo(t)%z=0的边界条件

y=1+t.*0;

return

function y=border_fune(t)%z=L的边界条件

y=t*.0;

return

运行情况：

按“run”运行时，弹出窗口

将图中相关数据更改为：

点击图框中的“OK”，在“command window”中输出结果为：