# How to plot a Phase-Portrait of three coupled differential equations?

I have the following problem: I would like to plot a StreamPlot (phase portrait) of three coupled, second-order, non-linear differential equations. These are:

m*l^2*θ1''[t] + m*x''[t]*l*Cos[θ1[t]] +
c1*θ1[t] + m*g*l*Sin[θ1[t]] == 0

m*l^2*θ2''[t] + m*x''[t]*l*Cos[θ2[t]] +
c1*θ2[t] + m*g*l*Sin[θ2[t]]==0

(M + 2 m)*x''[t] + c2*x'[t] +
k*x[t] + (m*l*θ1''[t]*Cos[θ1[t]] -
m*l*θ1'[t]^2*Sin[θ1[t]]) + (m*l*θ2''[t]*
Cos[θ2[t]] - m*l*θ2'[t]^2*Sin[θ2[t]])==0


where theta1(t), theta2(t) and x(t) are generalized coordinates and the rest are constants. Has anyone a hint how to plot a StreamPlot of e.g. theta1'(t) vs. theta1(t)?

Thank you very much. Best, Jonas.

• Try to solve your ode-system using ParametricNDSolve with variying initial conditions as parameter Commented Jul 21, 2022 at 8:54
• If you want to plot a phase plot, or any plot in general, you have to pass a number to all of the parameters (m, g, c1 etc.)
– tush
Commented Jul 21, 2022 at 16:16

Not sure what you want but here's a start: Can solve these numerically as IVPs if given values to all constants and starting conditions. Here's a set up with random values:

   (* set all constants to some vals *)
m = 1;
l = 1;
c1 = 1;
c2 = 1;
g = 1
M = 1
k = 1;
(* create an array of equations *)
theEqns = {m*l^2*\[Theta]1''[t] + m*x''[t]*l*Cos[\[Theta]1[t]] +
c1*\[Theta]1[t] + m*g*l*Sin[\[Theta]1[t]] == 0,
m*l^2*\[Theta]2''[t] + m*x''[t]*l*Cos[\[Theta]2[t]] +
c1*\[Theta]2[t] + m*g*l*Sin[\[Theta]2[t]] == 0,
(M + 2 m)*x''[t] + c2*x'[t] +
k*x[t] + (m*l*\[Theta]1''[t]*Cos[\[Theta]1[t]] -
m*l*\[Theta]1'[t]^2*Sin[\[Theta]1[t]]) + (m*l*\[Theta]2''[t]*
Cos[\[Theta]2[t]] - m*l*\[Theta]2'[t]^2*Sin[\[Theta]2[t]]) == 0}
(* define some initial conditions for the IVP *)
theInit = {\[Theta]1[0] == 0, \[Theta]1'[0] == 0.1, \[Theta]2[0] ==
0.5, \[Theta]2'[0] == 0.17, x[0] == 0, x'[0] == 0.2}


Now use NDSolveValue to solve the system for 0<=t<=25:

{theta1, theta2, theX} =
NDSolveValue[
Join[theEqns, theInit], {\[Theta]1, \[Theta]2, x}, {t, 0, 25}]


Now take the (numeric) derivative of theta1 and plot it vs. theta1[t]:

theta1D[t_] = D[theta1[t], t];
ParametricPlot[{theta1D[t], theta1[t]}, {t, 0, 25}]


• Hi Josh, thank you very much for your answer, it helped me! Best, Jonas. Commented Jul 23, 2022 at 2:54