How to graph a complex system of ODEs?

6 views (last 30 days)
Diego Torres-Siclait
Diego Torres-Siclait on 25 Apr 2020
Commented: Steven Lord on 26 Apr 2020
Any suggestions on how to graph this system?
(1) dU/dt = C1 * ( 1 - C2*U)
(2) dV/dt = C2 * dU/dt + C3 * (V - W)
(3) dW/dt = C4 * (V - W)
Where the Cs are all unique constants
I tried using ode45 to create arrays but am having trouble figuring out how to weave together several differential MATLAB functions that all depend on one another especiallly given how equation 2 depends on equation 1 and equation 2 and 3 are intertwined with shared variables.
This is a system of equations describing an isothermal kinetics PFR reactor w/ Heat Exchanger in terms of simplified mathematical variables (rather than the many constants and X, T, Ta, and V variables). Normally PolyMath is used but that program is very buggy.

Answers (1)

James Tursa
James Tursa on 25 Apr 2020
Edited: James Tursa on 25 Apr 2020
You have three 1st order equations, so you need a 3-element state vector for this. To keep the same nomenclature as the MATLAB docs, I will use y for this state vector. The definitions would be
y(1) = U
y(2) = V
y(3) = W
Then given a 3-element state vector y, write a derivative function for ydot. E.g.,
function ydot = myderivative(t,y,C1,C2,C3,C4)
ydot = zeros(size(y));
ydot(1) =
ydot(2) =
ydot(3) =
end
You write code for the three lines above. Wherever you see a U you would use y(1). Wherever you see a V you would use y(2). Wherever you would see a W you would use y(3). And dU/dt is simply ydot(1), be sure to calculate that first so you can use it in ydot(2). So write that code. Then your ode45( ) call would look something like:
[T,Y] = ode45(@(t,y)myderivative(t,y,C1,C2,C3,C4,tspan,y0);
  3 Comments
Steven Lord
Steven Lord on 26 Apr 2020
To make the code in myderivative look even more similar to the equations, you can define some intermediate variables.
function ydot = myderivative(t,y,C1,C2,C3,C4)
U = y(1);
V = y(2);
W = y(3);
dUdt = % write it terms of t, U, V, W, etc.
dVdt =
dWdt =
ydot = [dUdt; dVdt; dWdt];
end

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!