2
$\begingroup$

I need to force a whenevent action to be evaluated numerically all the times the event happens. The fact is that I use very complicated matrices and I need to invert them so the system requires a lot of time to work. Is there a way to make whenevent work only with numerical values in order to go faster?
I post a sort of toy model to represent what would I want to do:

Attributes[WhenEvent] = {};
ser[t_] = {x1[t], x2[t]};
Aw = {{-1, 2}, {0, -4}};
fun := Inverse[{{x1[t]*x2[t]^2,  37}, {x1[t]*Cos[x2[t]], 2} }];
eqinn = ser[0] == {1, 1};

eqdyn = D[ser[t], t] == Aw.ser[t];
eqinn = ser[0] == {1, 1};

ev = WhenEvent[
  Mod[t, 1], {ser[t] ->  fun.ser[t], "RestartIntegration"}]

sol = NDSolveValue[{eqdyn, eqinn, ev}, ser[t], {t, 0, 4}];
Plot[sol, {t, 0, 4}]

enter image description here

Now If I execute ev I obtain:

ev

(*WhenEvent[
 Mod[t, 1], {{x1[t], 
    x2[t]} -> {(2 x1[t])/(-37 Cos[x2[t]] x1[t] + 2 x1[t] x2[t]^2) - (
     37 x2[t])/(-37 Cos[x2[t]] x1[t] + 2 x1[t] x2[t]^2), -((
      Cos[x2[t]] x1[t]^2)/(-37 Cos[x2[t]] x1[t] + 
       2 x1[t] x2[t]^2)) + (
     x1[t] x2[t]^3)/(-37 Cos[x2[t]] x1[t] + 2 x1[t] x2[t]^2)}, 
  "RestartIntegration"}]*)

Just because it evaluates fun. In the real problem I can't do this way. A solution is to use the pattern test _?NumericQ, but I ask if is there another way because the function arguments would be too much in too much functions so I wouldn't write them every time as function arguments

So my goal is to have:

ev

(*WhenEvent[
 Mod[t, 1], {{x1[t], 
    x2[t]} -> fun.ser[t]}, 
  "RestartIntegration"}]*)

I would want that every time event happens It substitutes symbolic values with numerical values and then go on with matrix inversion

Improve this question
$\endgroup$
2
  • $\begingroup$ Why do you need to invert the matrices? I suggest that all you want is to solve linear systems. You rarely need to actually invert a matrix which takes a lot of time as you noticed. $\endgroup$
    – Somos
    Apr 1, 2021 at 20:37
  • $\begingroup$ @Somos In the real task I need to implement an hybrid kalman filter, and for the update step I need to invert a matrix. Thanks for the advice I'll go deep with your suggestion $\endgroup$
    – PeaceEverybody
    Apr 2, 2021 at 9:08

2 Answers 2

Reset to default
3
$\begingroup$

For this system we can make numerical code, but we need to test it with a large system

Attributes[WhenEvent] = {};
ser[t_] = {x1[t], x2[t]};
Aw = {{-1, 2}, {0, -4}};

func = Block[{x, y}, 
   With[{code = Inverse[{{x*y^2, 37}, {x*Cos[y], 2}}]}, 
    Compile[{{x, _Real}, {y, _Real}}, code, 
     CompilationTarget -> "C"]]];
funQ[x_?NumericQ, y_?NumericQ] := func[x, y];

eqinn = ser[0] == {1, 1};

eqdyn = D[ser[t], t] == Aw . ser[t];
eqinn = ser[0] == {1, 1};


ev = WhenEvent[
   Mod[t, 1], {ser[t] -> funQ[x1[t], x2[t]] . ser[t], 
    "RestartIntegration"}];


sol = NDSolveValue[{eqdyn, eqinn, ev}, ser[t], {t, 0, 4}];

Plot[sol, {t, 0, 4}, PlotRange -> All]

Figure 1

Improve this answer
$\endgroup$
7
  • $\begingroup$ I didn't know the use of compile. Thanks for the suggestion, I'll try it in my original code. $\endgroup$
    – PeaceEverybody
    Apr 2, 2021 at 9:10
  • $\begingroup$ @PeaceEverybody In a case of singular matrix when Inverse and LinearSolve not working we can use PseudoInverse with this code. $\endgroup$
    – Alex Trounev
    Apr 2, 2021 at 10:42
  • $\begingroup$ If I understand this right, this code computes the inverse symbolically and then compiles the result, where the code of @PeaceEverybody just uses the symbolic inverse. So the main difference is that this code compiles func -- is that right? (BTW, you can get rid of funQ if you use the Compile option RuntimeOptions -> {"EvaluateSymbolically" -> False}.) $\endgroup$
    – Michael E2
    Apr 2, 2021 at 13:38
  • $\begingroup$ @MichaelE2 Yes, your are right that with RuntimeOptions -> {"EvaluateSymbolically" -> False} we can avoid message from Compile. $\endgroup$
    – Alex Trounev
    Apr 2, 2021 at 14:38
  • $\begingroup$ @AlexTrounev in my original code I have all the variables time dependent, do you know if is there a way to force the numerical evaluation when the variable is function of x[t] and not x ? For example I got something like A which is a 5X5 matrix full of time dependent variables and I want it is evaluated only when it is numeric inside the ndsolve in which the time dependent variables are calculated $\endgroup$
    – PeaceEverybody
    Apr 6, 2021 at 9:58
3
$\begingroup$

I'd use LinearSolve instead of Inverse:

ser[t_] = {x1[t], x2[t]};
Aw = {{-1, 2}, {0, -4}};
mat := {{x1[t]*x2[t]^2, 37}, {x1[t]*Cos[x2[t]], 2}}; (*** No Inverse[] ***)
eqinn = ser[0] == {1, 1};

eqdyn = D[ser[t], t] == Aw . ser[t];
eqinn = ser[0] == {1, 1};

ev = WhenEvent[
      Mod[t, 1], {#1 -> LinearSolve[#2, #1], "RestartIntegration"}
      ] &[ser[t], mat];

sol = NDSolveValue[{eqdyn, eqinn, ev}, ser[t], {t, 0, 4}];
Plot[sol, {t, 0, 4}]

ev
(*
WhenEvent[
 Mod[t, 1], {{x1[t], x2[t]} -> 
   LinearSolve[{{x1[t] x2[t]^2, 37}, {Cos[x2[t]] x1[t], 2}}, {x1[t], 
     x2[t]}], "RestartIntegration"}]
*)
Improve this answer
$\endgroup$
3
  • $\begingroup$ really thanks for showing me the LinearSolve command, I'm gonna try to see if I'm do use it on my original code ! $\endgroup$
    – PeaceEverybody
    Apr 2, 2021 at 9:14
  • $\begingroup$ Nice code (+1). I think that your code could be effective in a case of large matrix. Only you need to replace LinearSolve with PseudoInverse in a case of singular matrix. $\endgroup$
    – Alex Trounev
    Apr 2, 2021 at 14:51
  • 1
    $\begingroup$ @AlexTrounev Thanks. I'd suggest LeastSquares, which is to PseudoInverse as LinearSolve is to Inverse. $\endgroup$
    – Michael E2
    Apr 2, 2021 at 15:03

Your Answer

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

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