Solving first order ODE with initial conditions and symbolic function

2 views (last 30 days)
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 Comments
Dyuman Joshi
Dyuman Joshi on 1 Dec 2023
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
Valerie on 1 Dec 2023
Above if the ODE I'm attempting to solve ^.
Output from Mathematica alongside code is:

Sign in to comment.

Accepted Answer

Torsten
Torsten on 1 Dec 2023
Edited: Torsten on 1 Dec 2023
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))
Tq = 3.3454e+02 + 1.2336e-71i
  4 Comments
Torsten
Torsten on 1 Dec 2023
Edited: Torsten on 1 Dec 2023
It asks to derive the temperature Tq at time tq = 2000 from the solution Tsol(2).

Sign in to comment.

More Answers (1)

Alan Stevens
Alan Stevens on 1 Dec 2023
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')

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!