# solving a coupled differential equations using ODE45

29 views (last 30 days)
Rithwik Kandukuri on 3 Dec 2020
Answered: Bjorn Gustavsson on 3 Dec 2020
I've a simple question about solving a differential equation using ode45. my only issue is that there is another differential equation in it and im not sure if what i am doing is wrong here is my code and im also attaching my differential functions.
.The initial condition for work done is w0=0.l i cant just write [th,W] = ode45(@work, range,w0)[gives me all 0s as output] because it is expressed in terms of function P(which also has its own initial condition P1) right ? i just need help in understanding how to write the function for work.
P.S. I am aware that the pressure and work functions are similar. i just seperated them for troubleshooting sakes.
PFA the equations
global theta r bore stroke P1 T1 y shaft rod epsilon V1 n a theta_s theta_d Qin range
[th,P] = ode45(@pressure,range,P1);
plot(th,P)
function [V,dV] = vol(theta)
global r ;
V = (1/r) + ((r-1)*(1-cos(theta)))/2*r ;
dV = ((r-1)/2*r )*sin(theta);
end
function [dxb] =burnfrac(theta)
global theta_s theta_d n a ;
xb = 1- exp(-a*(((theta-theta_s)*(1/theta_d))^n ));
if theta>theta_s && theta<(theta_s+theta_d)
dxb= n*(a/theta_d)*(1-xb)*((theta-theta_s)/theta_d)^(n-1);
else
dxb=0;
end
end
function [dP]= pressure(theta,P)
global Qin y ;
[V ,dV] = vol(theta);
dxb=burnfrac(theta);
dQ= Qin *dxb ;
dP= - y* P/V * dV + ((y-1)/V) * dQ;
dW= P*dV;
end
function [dW]= work(theta,P)
global Qin y ;
[V ,dV] = vol(theta)
dxb=burnfrac(theta);
dQ= Qin *dxb ;
dP= - y* P/V * dV + ((y-1)/V) * dQ
dW= P*dV
end

Bjorn Gustavsson on 3 Dec 2020
You have to put both ODEs together and integrate the pair at once. Let's look at an example equation:
That pair of equations can be simultaneosly integrated if we put them into one ode-function:
function dQdtdPdt = odePair(t,QP)
P = QP(2);
Q = QP(1); % for readability
dQdt = -P+f(t); % whatever you need for f...
dPdt = -P*Q + f;
dQdtdPdt = [dQdt;dPdt];
end
That function you now integrate with ode45 (for example):
Q0P0 = [12,32];
t_span = [0 34];
QP = ode45(@(t,QP) odePair(t,QP),t_span,Q0P0);
Which will give you both Q(t) and P(t).
HTH

R2020b

### Community Treasure Hunt

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

Start Hunting!