24小时热门版块排行榜

返回列表

zhaohuiyuan

金虫 (著名写手)

应助: 2 (幼儿园)
金币: 1401.3
散金: 3625
红花: 21
帖子: 1143
在线: 412.6小时
虫号: 2005943
注册: 2012-09-16
性别: MM
专业: 工程热物理相关交叉领域

[求助] 一维热传导偏微分方程的数值解法求助

各位虫友们，你们好！由于我对编程软件运用的不是很熟悉，希望得到各位的帮助，我编出来的程序怎么都运行不出来不知道是为什么，希望各位高手能指点一下，我要求解的方程是（见图片），如果哪位能帮上我，让我的程序能运行出结果来的话，我将用所有的金币来酬谢，若能给出有用意见的也行，我将酌情给金币，谢谢！
程序如下
function heat_conduction() %一维齐次热传导方程
options={'空间半径ro','空间点数N' ,'时间点数M','密度p','比热c','导热系数s',‘'对流传热系数h',};
topic='seting';
lines=1;
def={'0.032','100','1000','1714.8','1.48','0.535','61.27'};
b=inputdlg(options,topic,lines,def);%创建输入对话
ro=eval(b{1});%eval执行字符串
N=eval(b{2});
M=eval(b{3});
p=eval(b{4});
c=eval(b{5});
s=eval(b{6});
h=eval(b{7});
%***************************************************
Fo=s/(p*c);
or=ro/N;%空间步长
r=0

r:ro;
r=r';
ot=Fo*or^2/a;%时间步长
tm=M*ot;%热传导的总时间tm
t=0

t:tm;
t=t';
%计算初值和边值
T=zeros(N+1,M+1);
Ti=init_fun(r);%fun为目标函数的表达式字符串或MATLAB自定义函数的函数柄
To=border_funo(t);
Te=border_fune(t);
T(:,1)=Ti;%t=0时的初始温度
T(1,

=To;%r=0时的温度值
T(N+1,

=Te;%r=r0时的温度
%用差分法求出温度T与半径r、时间t的关系
for k=1:M
m=2;
while m<=N&&Fo<=(ro/(2*r+or))
T(m,k+1)=(Fo+(or*Fo)/ro)*(T(m+1,k)+(1-(or*Fo)/r+2*Fo)*T(m,k))+Fo*T(m-1,k);
m=m+1;
end;
end;
%设置立体网格
for i=1:M+1
X(:,i)=r;
end;
for j=1:N+1
Y(j,

=t;
end
mesh(X,Y,T);
view([1 -1 1]);
xlabel('r');
ylabel('t');
zlabel('T');
function y=init_fun(r)%初值条件t=0
y=14+t.*0
return

function y=border_funo(t)%r=0的边界条件
for k=1:M
m=1;
while Fo<=0.5
T(m+1,k+1)=(1-2Fo)*T(m+1,k)+2Fo*T(m,k);
end;
end;
return

function y=border_fune(t)%r=ro的边界条件
Tf=40+t.*0;
A=(h*ot)/(p*c*or);
for k=1:M
m=N;
while (1-A-Fo)>=0
T(m+1,k+1)=A*Tf+(1-A-Fo)*T(m+1,k))+Fo*T(m,k);
end;
end;
return