# Solving system of ODEs with ode45 and multiple time dependent variables

27 views (last 30 days)
Homero Vazquez on 25 Feb 2021
Commented: Homero Vazquez on 26 Feb 2021
I was given a system of 2 ODE's to solve for 2 different temperatures in the span of 24 hours. I was able to solve for the temperature for a constant Temperature Ta and constant I but in reality these values change for every hour of the day. I was given a table with the values of Ta and I for every hour t. I do not know how to implement the values into my function for every hour t. I am fairly new to matlab so I'm sure somebody out there can help. Here is my code so far. I will also provide a picture of the table I was given.
syms I h1 Tw Tg Mw Cw dTwdt Ta h2
eqn1 = 0.043 * I == h1 * (Tw-Tg) +(Mw*Cw*dTwdt)
DTwdtEqn = solve(eqn1,dTwdt)
eqn2 = 0.048 * I + h1*(Tw-Tg) == h2*(Tg-Ta)
TgEqn = solve(eqn2,Tg)
dTgdt = diff(TgEqn)
tRange = [0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23]
T0 = [293; 298];
[tSol, TSol] = ode45(@ODEfun,tRange,T0);
waterTemp = TSol(:,1);
glassTemp = TSol(:,2);
plot(tSol,waterTemp,"LineWidth",2)
hold on
plot(tSol,glassTemp,"LineWidth",2)
hold off
function dTdt = ODEfun(t,T)
Ta = 200;
I = 20;
Mw = 5;
Cw = 4180;
h1 = 30;
h2 = 100;
Tw = T(1);
Tg = T(2);
TwPrime = (I.*(43/100) +h1 *(Tg - Tw))/(Mw*Cw);
TgPrime = Tw/(h1+h2) - ( I.*(6/125) + Ta.* h2 + Tw *h1)/(h1+h2)^2;
dTdt = [TwPrime;TgPrime];
end

Bjorn Gustavsson on 25 Feb 2021
In your ODE-fun you can "simply" interpolate your time-dependent coefficients. Something like this:
function dTdt = ODEfun(t,T,t4Ta,Ta,t4I,I)
Ta = interp1(t4Ta,Ta,t,'pchip'); % Sensiblest spline-interpolation scheme, I think
I = interp1(t4I,I,t,'pchip');
Mw = 5;
Cw = 4180;
h1 = 30;
h2 = 100;
Tw = T(1);
Tg = T(2);
TwPrime = (I.*(43/100) +h1 *(Tg - Tw))/(Mw*Cw);
TgPrime = Tw/(h1+h2) - ( I.*(6/125) + Ta.* h2 + Tw *h1)/(h1+h2)^2;
dTdt = [TwPrime;TgPrime];
end
Then you send in the relevant time-parameter-pairings for t4Ta, Ta, t4I and I:
[tSol, TSol] = ode45(@(t,T) ODEfun(t,T,0:23,Ta,0:23,I),tRange,T0);
HTH
Homero Vazquez on 26 Feb 2021
Thanks !

R2020b

### Community Treasure Hunt

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

Start Hunting!

Translated by