29 views (last 30 days)

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

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

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

Start Hunting!
## 0 Comments

Sign in to comment.