7
$\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"}]
$\endgroup$
6
  • 1
    $\begingroup$ I wonder if it might work with the rhs expressed using Piecewise? $\endgroup$ Oct 8, 2018 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$ Oct 8, 2018 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, 2018 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$ Oct 8, 2018 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, 2018 at 15:47

2 Answers 2

5
$\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

$\endgroup$
1
  • $\begingroup$ this approach works bueno under version 10.2. nice idea with the quadratic form. thanks! $\endgroup$
    – freddy90
    Oct 9, 2018 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.

$\endgroup$
4
  • $\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, 2018 at 16:18
  • 2
    $\begingroup$ v11.3 -- did you try a fresh kernel? $\endgroup$
    – Chris K
    Oct 8, 2018 at 16:25
  • $\begingroup$ jep, tried that but unfortunately doesn't change anything… $\endgroup$
    – freddy90
    Oct 8, 2018 at 16:29
  • $\begingroup$ Yeah, my colleague's v10.3 has the same problem. Bug? $\endgroup$
    – Chris K
    Oct 8, 2018 at 17:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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