6
$\begingroup$

I'm curious whether it is possible to solve switched linear systems within the framework of NDSolve. For example a system of linear ode's like $$x'(t) = \left\{\begin{array}{ll} A_1 x(t),& \text{if} \,\, x_1x_2\leq 0 \\ A_2 x(t), & \text{if} \,\, x_1x_2>0 \end{array}\right.$$

where $A_1$ and $A_2$ are two constant matrices with appropriate size (namely $A_1,A_2 \in \mathbb{R}^{2\times 2}$) and $x(t) = \left(x_1(t),x_2(t)\right)^\top$.

I tried WhenEvent but received an error message saying

"Warning: the rule !(*SuperscriptBox[\"x\", \"[Prime]\", MultilineFunction->None][t] -> A1 . x[t]) will not directly set the state because the left-hand side is not a list of state variables."

Here is the code

A1 = {{0, -1}, {2, 0}};
A2 = {{0, -2}, {1, 0}};
x[t_] = {x1[t], x2[t]}

NDSolve[{x'[t] == A2.x[t], x1[0] == 6, x2[0] == 3, 
WhenEvent[x1[t] x2[t] <= 0, x'[t] -> A1.x[t]]}, {x1, x2}, {t, 0, 
100}, Method -> {"EquationSimplification" -> "Residual"}]
|improve this question
$\endgroup$
  • 1
    $\begingroup$ I wonder if it might work with the rhs expressed using Piecewise? $\endgroup$ – Daniel Lichtblau Oct 8 '18 at 14:50
  • $\begingroup$ Indeed, if memory serves, NDSolve[] will set up the WhenEvent[] objects on your behalf if you use Piecewise[]. Still, it is useful to know how to adapt WhenEvent[] in case the automatic method fails. $\endgroup$ – J. M. will be back soon♦ Oct 8 '18 at 15:01
  • 1
    $\begingroup$ Tried this: system = x'[t] == Piecewise[{{A1.x[t], x1[t] x2[t] < 0}, {A2.x[t], x1[t] x2[t] > 0}}]; NDSolve[{system, x1[t] == 3, x2[0] == 2}, {x1, x2}, {t, 0, 10}]. I got the error: "Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions" $\endgroup$ – freddy90 Oct 8 '18 at 15:02
  • 1
    $\begingroup$ Why are the initial conditions for both x1 and x2 scalars and not vectors? Try using Indexed[] if you want to refer to a vector-valued function componentwise, just like in your inequality conditions. $\endgroup$ – J. M. will be back soon♦ Oct 8 '18 at 15:05
  • 1
    $\begingroup$ It doesn't seem that defining x[t_] = {x1[t], x2[t]} is enough for Mathematica to know the relationship between x[t] and {x1[t],x2[t]}. $\endgroup$ – Chris K Oct 8 '18 at 15:47
4
$\begingroup$

You can use Piecewise in the vector form of the ODE, the only tricky part is how to create the condition. Here are two possibilities:

pm1 = {{0, 1}, {1, 0}};

sol1 = NDSolveValue[
    {
    x'[t] == Piecewise[{{A2.x[t], x[t].pm1.x[t]>0}}, A1.x[t]],
    x[0] == {6, 3}
    },
    x,
    {t, 0, 100}
];

sol2 = NDSolveValue[
    {
    x'[t] == Piecewise[{{A2.x[t], Indexed[x[t], 1] Indexed[x[t], 2] > 0}}, A1.x[t]],
    x[0] == {6, 3}
    },
    x,
    {t, 0, 100}
];

Visualizations:

Plot[
    {Indexed[sol1[t],1], Indexed[sol1[t], 2]},
    {t,0,100},
    PlotRange->All
]

enter image description here

Plot[
    {Indexed[sol2[t],1], Indexed[sol2[t], 2]},
    {t,0,100},
    PlotRange->All
]

enter image description here

|improve this answer
$\endgroup$
  • $\begingroup$ this approach works bueno under version 10.2. nice idea with the quadratic form. thanks! $\endgroup$ – freddy90 Oct 9 '18 at 15:12
6
$\begingroup$

Like the error message says, I believe WhenEvent needs to change a state variable, not its (highest) derivative. Here's an approach that sets A as a DiscreteVariable that can be changed when needed.

listProduct[x_List] := Times @@ x;

A1 = {{0, -1}, {2, 0}};
A2 = {{0, -2}, {1, 0}};

sol = NDSolve[{x'[t] == A[t].x[t], x[0] == {6, 3}, A[0] == A2,
  WhenEvent[listProduct[x[t]] <= 0, A[t] -> A1],
  WhenEvent[listProduct[x[t]] > 0, A[t] -> A2]}, {x, A}, {t, 0, 100},
  DiscreteVariables -> {A}][[1]];

Plot[Sign[listProduct[x[t] /. sol]], {t, 0, 100}]

Mathematica graphics

listProduct is by rm -rf from this answer.

|improve this answer
$\endgroup$
  • $\begingroup$ What version of mathematica do you use? I tried your code with version 10.2 and I can only integrate till t = 4 without waiting forever $\endgroup$ – freddy90 Oct 8 '18 at 16:18
  • 2
    $\begingroup$ v11.3 -- did you try a fresh kernel? $\endgroup$ – Chris K Oct 8 '18 at 16:25
  • $\begingroup$ jep, tried that but unfortunately doesn't change anything… $\endgroup$ – freddy90 Oct 8 '18 at 16:29
  • $\begingroup$ Yeah, my colleague's v10.3 has the same problem. Bug? $\endgroup$ – Chris K Oct 8 '18 at 17:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.