Solving first order ODE with initial conditions and symbolic function
显示 更早的评论
The code returns a solution involving a complex number. I know this is not correct because I have solved it in Mathematica. Is there a way to solve it in MATLAB? Below is my code with all variables defined:
% input parameters
Tinf=70+273.15;
Ti=20+273.15;
d=15e-2;
r=d/2;
cdepth=10/1100;
Tf=50+273.15;
cp=4183;
rho=994;
g=9.81;
mew=0.007196;
k=0.6107;
Pr=4.929
Beta=0.00347;
L=10e-2;
% define equations
kv=mew/rho;
Asc=pi*r^2;
V=pi*r^2*L;
m=rho*cp*V;
% Calculate the Ray Number
Gr =@(T) (g*Beta*(Tinf-T)*L^(3))/(kv^(2))
Ray =@(T) Gr(T)*Pr
Nu =@(T) 0.15*(Ray(T)^(1/3))
h =@(T) (Nu(T)*k)/(L)
syms T(t) ;
ode = m*diff(T) == Asc*h(T)*(Tinf-T)
cond = T(0) == Ti;
TSol(t) = dsolve(ode,cond)
disp(TSol)
2 个评论
Dyuman Joshi
2023-12-1
Please share the mathematical definition of ODE that you are trying to solve.
Also, please share the output from Mathematica, along with the code used there.
Valerie
2023-12-1
Above if the ODE I'm attempting to solve ^.
Output from Mathematica alongside code is:
采纳的回答
You solved it in your code. The second solution out of the three MATLAB returned is the "correct" one giving real-valued temperatures. You can ignore the complex component of size 1e-71. The three solutions result from the T^1/3 term in the Nusselt number.
% input parameters
Tinf=70+273.15;
Ti=20+273.15;
d=15e-2;
r=d/2;
cdepth=10/1100;
Tf=50+273.15;
cp=4183;
rho=994;
g=9.81;
mew=0.007196;
k=0.6107;
Pr=4.929;
Beta=0.00347;
L=10e-2;
% define equations
kv=mew/rho;
Asc=pi*r^2;
V=pi*r^2*L;
m=rho*cp*V;
syms T(t) t
% Calculate the Ray Number
Gr = g*Beta*(Tinf-T)*L^3/kv^2;
Ray = Gr*Pr;
Nu = 0.15*Ray^(1/3);
h = Nu*k/L;
ode = m*diff(T) == Asc*h*(Tinf-T);
Tsol = dsolve(ode,T(0)==Ti);
fplot(real(Tsol(2)),[0 10000])
tq = 2000;
Tq = double(subs(Tsol(2),t,tq))
更多回答（1 个）
You can do it numerically as follows:
% input parameters
Tinf=70+273.15;
Ti=20+273.15;
d=15e-2;
r=d/2;
cdepth=10/1100;
Tf=50+273.15;
cp=4183;
rho=994;
g=9.81;
mew=0.007196;
k=0.6107;
Pr=4.929;
Beta=0.00347;
L=10e-2;
% define equations
kv=mew/rho;
Asc=pi*r^2;
V=pi*r^2*L;
m=rho*cp*V;
% Calculate the Ray Number
Gr =@(T) g*Beta*(Tinf-T)*L^3/kv^2;
Ray =@(T) Gr(T)*Pr;
Nu =@(T) 0.15*Ray(T)^(1/3);
h =@(T) Nu(T)*k/L;
dTdt = @(t,T)Asc/m*h(T)*(Tinf-T);
tend = 10^4;
tspan = [0 tend];
[t,Tsol] = ode45(dTdt, tspan, Ti);
plot(t,Tsol,[0 tend],[Tinf Tinf],'--'), grid
xlabel('t'), ylabel('T')
类别
在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Translated by 
