# Solving first order ODE with initial conditions and symbolic function

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

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
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).
Valerie on 2 Dec 2023
Thank you so much!

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')
Valerie on 1 Dec 2023
Is there a way to do this symbolically since I'm trying to solve for a particular time?