1
$\begingroup$

I am solving a system of ODE which contains a discontinuous ode (the equation v[t]== ...in the following code and it means that my v[t] jumps between vmax (i.e.,2.25) and 0).

After a long time of trial and error (try to solve the odes system by using If, PieceWise, UnitStep, etc), I found that WhenEvent successfully gives me a reasonable result (sol1 in my code).

Then I tried to solve my system of ode by using smooth function of (Tanh) to approximate the discontinuous ode, and then I manually converted the algebraic equation to an ode by taking derivative with respect to time as in this post. This time I solved the odes much faster, but I got a very different and unreasonable result (sol2 in my code).

My code is shown below.

Why are these results (i.e., sol1 by using event and sol2 by converting to odes ) so different? What's wrong with my sol2? Is it beacuse of my converting the algebraic equation of v[t] to a ode? Or Tanh function not approroate here?

Thanks.


My code:

Remove["Global`*"] // Quiet;

tend = 2;
tdat = Range[0, tend, 0.01];
gdat =(*RandomVariate[NormalDistribution[0,0.5],Length[tdat]]*)
  2 Exp[-0.5 tdat] Sin[100 tdat]^2 Cos[50 tdat]^3;
tgdat = {tdat, gdat}\[Transpose];
xg = Interpolation[tgdat, InterpolationOrder -> 1];
(*Plot[xg[t],{t,0,tend},PlotRange\[Rule]All];*)

m = {{98.3, 0, 0}, {0, 98.3, 0}, {0, 0, 98.3}};
c = {{175, -50, 0}, {-50, 100, -50}, {0, -50, 50}};
k = {{12, -6.84, 0}, {-6.84, 13.7, -6.84}, {0, -6.84, 6.84}} 10^5;
x = {x1[t], x2[t], x3[t]};
dx = {x1'[t], x2'[t], x3'[t]};
ddx = {x1''[t], x2''[t], x3''[t]};

Α = 301;
n = 2;
η = 190;

vmax = 2.25;
αs[u_] := 140 10^2 + 695 10^2 u;
c1[u_] := 283 10^2 + 2.95 10^2 u;
c0[u_] := 21 10^2 + 3.5 10^2 u;
f[u_, dy_, x1_] := c1[u] dy + 500 (x1 - 0);
fc[x1_, x2_, x3_, v1_, v2_, v3_] := 
  4 10^6 x1 - 3.4 10^6 x2 + 752341 x3 + 37284 v1 + 26262 v2 + 7604 v3;

sol1 = NDSolve[{
      Flatten[m.ddx + c.dx + k.x] == 
       Flatten[-m.{1, 1, 1} xg[t] + {-1, 0, 0} f[u[t], y'[t], 
           x1[t]]],
      z'[t] == -50 10^4 Sqrt[(x1'[t] - y'[t])^2] z[t] Sqrt[z[t]^2] - 
        363 10^4 (x1'[t] - y'[t]) z[
          t]^2 + Α (x1'[t] - y'[t]),
      y'[t] == 
       1/(c0[u[t]] + 
         c1[u[t]]) (αs[u[t]] z[t] + c0[u[t]] D[x1[t], t] + 
          4690 (x1[t] - y[t])),
      WhenEvent[(fc[x1[t], x2[t], x3[t], x1'[t], x2'[t], x3'[t]] - 
           f[u[t], y'[t], x1[t]]) f[u[t], y'[t], x1[t]] < 0, 
       v[t] -> 0],
      WhenEvent[(fc[x1[t], x2[t], x3[t], x1'[t], x2'[t], x3'[t]] - 
           f[u[t], y'[t], x1[t]]) f[u[t], y'[t], x1[t]] > 0, 
       v[t] -> 2.25],
      u'[t] == -η (u[t] - v[t]),
      v[0] == 0, u[0] == 0, z[0] == 0, y[0] == 0, 
      x1[0] == x1'[0] == 0, x2[0] == x2'[0] == 0, 
      x3[0] == x3'[0] == 0} // Flatten,
    {y, z, v, u, x1, x2, x3}, {t, 0, tend}
    , Method -> {"EquationSimplification" -> "Residual"}
    , MaxSteps -> Infinity
    , DiscreteVariables -> v[t]] // Flatten;
sol2 = NDSolve[{
      Flatten[m.ddx + c.dx + k.x] == 
       Flatten[-m.{1, 1, 1} xg[t] + {-1, 0, 0} f[u[t], y'[t], 
           x1[t]]],
      z'[t] == -50 10^4 Sqrt[(x1'[t] - y'[t])^2] z[t] Sqrt[z[t]^2] - 
        363 10^4 (x1'[t] - y'[t]) z[
          t]^2 + Α (x1'[t] - y'[t]),
      y'[t] == 
       1/(c0[u[t]] + 
         c1[u[t]]) (αs[u[t]] z[t] + c0[u[t]] D[x1[t], t] + 
          4690 (x1[t] - y[t])),
      v'[t] == 
       D[vmax (1 + 
          Tanh[2000 (fc[x1[t], x2[t], x3[t], x1'[t], x2'[t], x3'[t]] -
               f[u[t], y'[t], x1[t]]) f[u[t], y'[t], x1[t]]])/2, t],
      (*v[t]=vmax Simplify`PWToUnitStep@Piecewise[{{0,(fc[x1[t],x2[t],
      x3[t],x1'[t],x2'[t],x3'[t]]-f[u[t],y'[t],x1[t]])f[u[t],y'[t],x1[
      t]]≤ 0}},vmax],*)

      u'[t] == -η (u[t] - v[t]),
      v[0] == 0, u[0] == 0, z[0] == 0, y[0] == 0, 
      x1[0] == x1'[0] == 0, x2[0] == x2'[0] == 0, 
      x3[0] == x3'[0] == 0} // Flatten,
    {y, z, v, u, x1, x2, x3}, {t, 0, tend}
    , MaxSteps -> Infinity] // Flatten;
Plot[u[t] /. {sol1, sol2} // Evaluate, {t, 0, tend}, PlotRange -> All,
  PlotTheme -> "Web", PlotLegends -> {"using event", "using tanh"}]
Plot[v[t] /. {sol1, sol2} // Evaluate, {t, 0, tend}, PlotRange -> All,
  PlotTheme -> "Web", PlotLegends -> {"using event", "using tanh"}]
Plot[y[t] /. {sol1, sol2} // Evaluate, {t, 0, tend}, PlotRange -> All,
  PlotTheme -> "Web", PlotLegends -> {"using event", "using tanh"}]

enter image description here

enter image description here

$\endgroup$
  • 1
    $\begingroup$ Change your second WhenEvent v[t]->2.25 to v[t]->1 ! Now the solutions sol1 an sol2 are the same! $\endgroup$ – Ulrich Neumann Nov 9 '19 at 14:05
  • $\begingroup$ Thanks. I tried your suggestion, but sol1 is still apparently different from sol1 (for example u[t] of sol1 and sol2). BTW, why did you suggest that changing v[t]->2.25 to v[t]->1 ? my vmax is 2.25, not 1, and it jumps between 0 and vmax(2.25). $\endgroup$ – xinxin guo Nov 9 '19 at 14:12
  • $\begingroup$ Sorry I didn't recognize vmax in sol2. Setting v[t]->1gives a good match of sol1 and sol2 in your second last plot (don't know why) $\endgroup$ – Ulrich Neumann Nov 9 '19 at 14:17
  • $\begingroup$ Thanks. Acctually the most important result to me is u[t]. $\endgroup$ – xinxin guo Nov 9 '19 at 14:22
  • $\begingroup$ You might want to plot fc[x1[t], x2[t], x3[t], x1'[t], x2'[t], x3'[t]] - f[u[t], y'[t], x1[t]]) f[u[t], y'[t], x1[t]] versus v[t] to see if every event was captured by NDSolve. $\endgroup$ – Michael E2 Nov 9 '19 at 17:33
3
$\begingroup$

These are different systems of equations. You can coordinate them by removing equation v'[t]== in the second system and replacing v[t] to If[t < 0.05050561921089747, 0, vmax (1 + Tanh[.01 (fc[x1[t], x2[t], x3[t], x1'[t], x2'[t], x3'[t]] - f[u[t], y'[t], x1[t]]) f[u[t], y'[t], x1[t]]])/2]. Then the solutions coincide up to t = 1.37, but then diverge.

tend = 2;
tdat = Range[0, tend, 0.01];
gdat =(*RandomVariate[NormalDistribution[0,0.5],Length[tdat]]*)
  2 Exp[-0.5 tdat] Sin[100 tdat]^2 Cos[50 tdat]^3;
tgdat = {tdat, gdat}\[Transpose];
xg = Interpolation[tgdat, InterpolationOrder -> 1];
(*Plot[xg[t],{t,0,tend},PlotRange\[Rule]All];*)

m = {{98.3, 0, 0}, {0, 98.3, 0}, {0, 0, 98.3}};
c = {{175, -50, 0}, {-50, 100, -50}, {0, -50, 50}};
k = {{12, -6.84, 0}, {-6.84, 13.7, -6.84}, {0, -6.84, 6.84}} 10^5;
x = {x1[t], x2[t], x3[t]};
dx = {x1'[t], x2'[t], x3'[t]};
ddx = {x1''[t], x2''[t], x3''[t]};

\[CapitalAlpha] = 301;
n = 2;
\[Eta] = 190;

vmax = 2.25;
\[Alpha]s[u_] := 140 10^2 + 695 10^2 u;
c1[u_] := 283 10^2 + 2.95 10^2 u;
c0[u_] := 21 10^2 + 3.5 10^2 u;
f[u_, dy_, x1_] := c1[u] dy + 500 (x1 - 0);
fc[x1_, x2_, x3_, v1_, v2_, v3_] := 
  4 10^6 x1 - 3.4 10^6 x2 + 752341 x3 + 37284 v1 + 26262 v2 + 7604 v3;

sol1 = NDSolve[{Flatten[m.ddx + c.dx + k.x] == 
       Flatten[-m.{1, 1, 1} xg[t] + {-1, 0, 0} f[u[t], y'[t], x1[t]]],
       z'[t] == -50 10^4 Sqrt[(x1'[t] - y'[t])^2] z[t] Sqrt[z[t]^2] - 
        363 10^4 (x1'[t] - y'[t]) z[t]^2 + \[CapitalAlpha] (x1'[t] - 
           y'[t]), 
      y'[t] == 
       1/(c0[u[t]] + c1[u[t]]) (\[Alpha]s[u[t]] z[t] + 
          c0[u[t]] D[x1[t], t] + 4690 (x1[t] - y[t])), 
      WhenEvent[(fc[x1[t], x2[t], x3[t], x1'[t], x2'[t], x3'[t]] - 
           f[u[t], y'[t], x1[t]]) f[u[t], y'[t], x1[t]] < 0, 
       v[t] -> 0], 
      WhenEvent[(fc[x1[t], x2[t], x3[t], x1'[t], x2'[t], x3'[t]] - 
           f[u[t], y'[t], x1[t]]) f[u[t], y'[t], x1[t]] > 
        0, v[t] -> vmax], u'[t] == -\[Eta] (u[t] - v[t]), 
      v[0] == 0, u[0] == 0, z[0] == 0, y[0] == 0, 
      x1[0] == x1'[0] == 0, x2[0] == x2'[0] == 0, 
      x3[0] == x3'[0] == 0} // Flatten, {y, z, v, u, x1, x2, x3}, {t, 
     0, tend}, Method -> {"EquationSimplification" -> "Residual"}, 
    MaxSteps -> Infinity, DiscreteVariables -> v[t]] // Flatten;
sol2 = NDSolve[{Flatten[m.ddx + c.dx + k.x] == 
       Flatten[-m.{1, 1, 1} xg[t] + {-1, 0, 0} f[u[t], y'[t], x1[t]]],
       z'[t] == -50 10^4 Sqrt[(x1'[t] - y'[t])^2] z[t] Sqrt[z[t]^2] - 
        363 10^4 (x1'[t] - y'[t]) z[t]^2 + \[CapitalAlpha] (x1'[t] - 
           y'[t]), 
      y'[t] == 
       1/(c0[u[t]] + c1[u[t]]) (\[Alpha]s[u[t]] z[t] + 
          c0[u[t]] D[x1[t], t] + 4690 (x1[t] - y[t])), u'[t] == -\[Eta] (u[t] - 
          If[t < 0.05050561921089747, 0, 
           vmax (1 + 
               Tanh[.01 (fc[x1[t], x2[t], x3[t], x1'[t], x2'[t], 
                    x3'[t]] - f[u[t], y'[t], x1[t]]) f[u[t], y'[t], 
                  x1[t]]])/2]), u[0] == 0, z[0] == 0, y[0] == 0, 
      x1[0] == x1'[0] == 0, x2[0] == x2'[0] == 0, 
      x3[0] == x3'[0] == 0} // Flatten, {y, z, u, x1, x2, x3}, {t, 0, 
     tend}, MaxSteps -> Infinity, 
    Method -> {"EquationSimplification" -> "Residual"}] // Flatten;
Plot[u[t] /. {sol1, sol2} // Evaluate, {t, 0, tend}, PlotRange -> All,
  PlotTheme -> "Web", PlotLegends -> {"using event", "using tanh"}]

Plot[y[t] /. {sol1, sol2} // Evaluate, {t, 0, tend}, PlotRange -> All,
  PlotTheme -> "Web", PlotLegends -> {"using event", "using tanh"}]

Figure 1

We can achieve a complete match if we put in the first system

WhenEvent[(fc[x1[t], x2[t], x3[t], x1'[t], x2'[t], x3'[t]] - 
     f[u[t], y'[t], x1[t]]) f[u[t], y'[t], x1[t]] < 10^-3, 
 v[t] -> 0], WhenEvent[(fc[x1[t], x2[t], x3[t], x1'[t], x2'[t], 
      x3'[t]] - f[u[t], y'[t], x1[t]]) f[u[t], y'[t], x1[t]] > -1, v[t] -> vmax]

Figure 2

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Did you mean "but then diverge"? It looks like an event was missed around t = 1.37 and a couple of times after that. Could be the reason for the divergence. $\endgroup$ – Michael E2 Nov 9 '19 at 17:30
  • $\begingroup$ Thanks. I found that If[t < 10^-5, .... also produced almost the same result as yours, so it seems that the two systems only different in the initial conditions. I think this more. $\endgroup$ – xinxin guo Nov 9 '19 at 18:22
  • $\begingroup$ @MichaelE2 You're right, I corrected for diverge. Obviously the difference in the missed event. Found that this is due to an error in calculation fc. $\endgroup$ – Alex Trounev Nov 9 '19 at 19:11
  • $\begingroup$ @xinxinguo See update to my answer. $\endgroup$ – Alex Trounev Nov 9 '19 at 19:12
  • $\begingroup$ @ Alex Trounev, thanks, your answer is very helpful and makes me find many more interesting aspects of NDSolve. I really appreciate your time helping many Mathematica fans. Have a nice day :) $\endgroup$ – xinxin guo Nov 10 '19 at 11:37

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.