How to solve a systems of ODE and Algebraic Equations

21 views (last 30 days)
Telema Harry on 28 Jan 2021
Commented: Alex Sha on 29 Jan 2021
I have a system of 3 nonlinear ODE and 2 nonlinear algebraic equations.
Please how can I solve these systems of equation.
ODE 45 can easily solve the ODE part. However, I don't know how to combine the solution from ODE45 and the algebraic equations.
Thank you.
jessupj on 28 Jan 2021
it sounds like what you're after is "how to solve a DAE" if the algebraic eqations constrain the solutions of the ODE part https://www.mathworks.com/help/matlab/math/solve-differential-algebraic-equations-daes.html
otherwise, if the algebraic equations aren't constraints (ie. they determine diagnostic variables), you probably want to solve the ODE and then solve the algebraic equations 'offline' using e.g. fsolve

jessupj on 28 Jan 2021
Edited: jessupj on 28 Jan 2021
looks like you've got a non-autonomous DAE.
with u=x(4) and y = x(5), you'd have:
dx(1) = -wh.*x(1) + wh.* x(5)
dx(2) = -wl.*x(2) + A.*sin(w.*t).* wl.*(x(5) - x(1))
dx(3) = K.*x(2)
0 = x(3) + A.* sin(w.*t) - x(4)
0 = 25 - (5 - x(4) ).^2 - x(5) % = 25 - (25 -10*x4 + x4^2) -x5 = x4*(10 -x4)-x5
and check this old post:
Telema Harry on 29 Jan 2021
Thank you for your help @jessupj and @James Tursa.
I have been able to run the scripts but there is something I am not doing correctly.
I want to be able to plot my x(1) , x(2), x(3), x(4),, x(5) againt time. But I have not been able to figure out how to do that.
% x(1) = eta (n)
% x(2) = psi (E) - that is the second differential equation
% x(3) = Uhat
% x(4) = u
% x(5) = J
tspan = 0:4;
iCon = [0; 0; 0 ; 0; 0];
M = [1 0 0 0 0; 0 1 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0];
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6 1e-6 1e-6]);
[t,y] = ode15s(@ESCdae,tspan,iCon)
plot (t,y)
function out = ESCdae(t,x)
f = 10;
A = 0.2;
w = 2*pi*f;
wh = 0.8;
wl = 0.2;
K = 5;
out = [-wh.*x(1) + wh.* x(5)
-wl.*x(2) + A.*sin(w.*t).* wl.*(x(5) - x(1))
K.*x(2)
x(3) + A.* sin(w.* t) - x(4)
x(4).*(10 - x(4)) - x(5)
]
end

Telema Harry on 28 Jan 2021
Edited: Telema Harry on 28 Jan 2021
Thank you James for the feedback. See the sample equations.
f = 10;
A = 0.2;
w = 2*pi*f;
wh = 0.8;
wl = 0.2;
y = 25 - (5 - u).^2
dx1dt = -wh.*x(1) + wh.* y
dx2dt = -wl.*x(2) + A.*sin(w.*t).* wl.*(y - x(1))
dx3dt = K.*x(2) % the solution to this equation gives x(3)
u = x(3) + A.* sin(w.*t)
Alex Sha on 29 Jan 2021
Hi, since: u = x(3) + A.* sin(w.*t) and y = 25 - (5 - u).^2, so y = 25 - (5 - ( x(3) + A.* sin(w.*t))).^2, substitute y into dx1dt and dx2dt, then pure ODE functions will be formed.

R2020a

Community Treasure Hunt

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

Start Hunting!