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)

[Dyuman Joshi]

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]

Above if the ODE I'm attempting to solve ^.

Output from Mathematica alongside code is:

Answer 1

[Torsten]

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

[Valerie]

Thank you so much!

[Valerie]

Is this part solving for a time? What is this code doing?

tq = 2000;
Tq = double(subs(Tsol(2),t,tq))

[Torsten]

It asks to derive the temperature Tq at time tq = 2000 from the solution Tsol(2).

[Valerie]

Thank you so much!

Answer 2

[Alan Stevens]

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]

Is there a way to do this symbolically since I'm trying to solve for a particular time?