MATLAB: How to solve a 3 ODEs with 1 vector input and constant parameters with 3 initial conditions parameters

differential equationsmultiple odeode45

Following this, an example of 3 ODEs matrix that I need to solve:
function Tdot=tempdyn(t,x)
global Tamb C1 C3 C4 T1 T3 T4 Wd Ws Wc % Only Wc is a vector input, other global variables are constants
Te=x(1); Ts=x(2); Tc=x(3);
Tedot=Te*(-1)*(1/(C4*T3)+1/(C4*T4))+Ts/(C4*T3)+(Tamb+273.15)/(C4*T4);
Tsdot=Te/(C3*T3)+Ts*(-1)*(1/(C3*T1)+1/(C3*T3))+Tc/(C3*T1)+(Wd/2+Ws)/(C3);
Tcdot=Ts/(C1*T1)-Tc/(C1*T1)+(Wd/2+Wc)./(C1);
Tdot=[Tedot;Tsdot;Tcdot];
Here, it is the main script to call the ode function. Only Wc is a vector input of 8760 and other parameters are constants and calculated previously. I use global to call each parameter
Wc=3*Rac*Ic.^2; % conductor losses
Ws=0; % screen losses assuming 0
% Initial conditions in Kelvin degrees
Te0=0+273.15;
Ts0=0+273.15;
Tc0=0+273.15;
tspan=[0 8760];
temp0=[Te0;Ts0;Tc0];
% Solve ODE using ODE numerical solver function ode45
[T,X]=ode45('tempdyn',tspan,temp0)
% Plotting
figure()
plot(T,X(:,3),'r-.')
xlabel('Time [hours]')
ylabel('Temperature [\circK]')
I want to solve this ode function formed by a 3-equation differential system. Any help would be very appreciated and thank you very much for your time consideration. I look forward to hearing back from you. Do not hesitate to contact me if you require further information.
Kind regards,

Best Answer

Wc=3*Rac*Ic.^2; % conductor losses
Ws=0; % screen losses assuming 0
% Initial conditions in Kelvin degrees
Te0=0+273.15;
Ts0=0+273.15;
Tc0=0+273.15;
tspan=[1 8760];
temp0=[Te0;Ts0;Tc0];
% Solve ODE using ODE numerical solver function ode45
[T,X]=ode45(@tempdyn,tspan,temp0);
% Plotting
figure();
plot(T,X(:,3),'r-.');
xlabel('Time [hours]');
ylabel('Temperature [\circK]'); 
function Tdot=tempdyn(t,x)
global Tamb C1 C3 C4 T1 T3 T4 Wd Ws Wc % Only Wc is a vector input, other global variables are constants
T=1:8760;
Wc_actual=interp1(T,Wc,t);
Te=x(1); 
Ts=x(2); 
Tc=x(3);
Tedot=Te*(-1)*(1/(C4*T3)+1/(C4*T4))+Ts/(C4*T3)+(Tamb+273.15)/(C4*T4);
Tsdot=Te/(C3*T3)+Ts*(-1)*(1/(C3*T1)+1/(C3*T3))+Tc/(C3*T1)+(Wd/2+Ws)/(C3);
Tcdot=Ts/(C1*T1)-Tc/(C1*T1)+(Wd/2+Wc_actual)/(C1);
Tdot=[Tedot;Tsdot;Tcdot];
Best wishes
Torsten.
Related Question