solving a coupled differential equations using ODE45
228 views (last 30 days)
Show older comments
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.
Thanks in advance.
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
0 Comments
Accepted Answer
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
0 Comments
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!