Help using WhenEvent with NDSolve to avoid complex solutions

Questions
,
1

Hello there, I would appreciate your help with me on this. I am modeling a second order differential equation ( mass and spring ). The spring is nonlinear (k x[t]^1.5). NDSolve solves the system but when encounters complex solutions it stoppes which makes sense to me. I want to make a constraint or a condition by using WhenEvent function inside the NDSolve to set x[t] to zero everytime the solution becomes complex or to just take the positive real part of it and continue solving for the specified period of time of 10 seconds. I’ve trying with no success. Please see the code in this message and I am also attaching the file for you. I would highly appreciate your help. Thank you

values = {m1 -> 1, m2 -> 2, m3 -> 3, E1 -> 210, E2 -> 210, 
   E3 -> 210, \[Nu]1 -> 1/3, \[Nu]2 ->  1/3, \[Nu]3 ->  1/3, 
   R1 -> 349/10000, R2 -> 349/10000, R3 -> 349/10000};
IC = {x1[0] == 0, x1'[0] == 1};
k1 = (1 - \[Nu]1^2)/(\[Pi] E1) /. values;
k2 = (1 - \[Nu]2^2)/(\[Pi] E2) /. values;
K12 = 4/(3 \[Pi] (k1 + k2)) Sqrt[(R1 R2)/(R1 + R2)] /. values;
\[Delta]12 = x1[t] + R1 /. values;
F12 = K12 \[Delta]12^1.5;
EQ1 = m1 x1''[t] +  K12 (x1[t])^(3/2) == 0 /. values;
eqn = NDSolve[{EQ1, IC} /. values, {x1[t],x1'[t]
}, {t, 0, 10}, Method -> "StiffnessSwitching"]
Plot[Evaluate[x1[t] /. eqn], {t, 0, 10}, GridLines -> Automatic, 
 PlotRange -> All]
Plot[Evaluate[x1'[t] /. eqn], {t, 0, 10}, GridLines -> Automatic, 
 PlotRange -> All]
2

The equation from your notebook is:

EQ1 = m1 x1''[t] + (x1[t])^(3/2) == 0 /. values

You begin with x1’[t] and x1[t] positive, but the force term drives x1’[t] negative and eventually x1[t] goes negative. This is when (x1[t])^(3/2) becomes complex. If this is a projectile, it just hit the ground and you could use “StopIntegration” in WhenEvent.

If it is really a spring which needs a real restoring force when x1 < 0, then the mathematics needs a different formulation. For example:

EQ1 = m1 x1''[t] + Sign[x1[t]] Abs[x1[t]]^(3/2) == 0 /. values

which calculates a real spring force, but also preserves the sign of the force.

Kind regards,
David

3

Dear David, thank you for your help. I understand your explanation of the problem. It made it a lot easier for me to understand what the NDSolve does when it encounters negative solutions, which results in complex solutions. Would you please tell me what are you telling NDSolve to do by your suggestion ( Sign[x1[t]] Abs[x1[t]]^(3/2) ). I want to understand what does it do when it encounters negative solutuions.
Thank you again for your help David
Regards, Khalid

4

Hi Khalid,

I rewrite your original equation a bit:

In[11]:= eq1 = m1 x''[t] == - (x[t])^((3/2))

We see the right hand term is a restoring force. When x[t] is positive, it applies a negative force to send x in the negative direction. If it’s a real spring, then when x[t] goes negative, it should be applying a positive force. But it doesn’t. Instead it produces a complex number.

So I changed the physics:

eq2 = x''[t] == -Sign[x[t]] Abs[x[t]]^(3/2)

Now the term on the right hand side works for both positive and negative x[t]. It has the same nonlinear relation on each side of zero, with the proper change of sign. Look up the functions Sign and Abs. Here we plot the restoring force term.

Plot[-Sign[x[t]] Abs[x[t]]^(3/2), {x[t], -1, 1}, 
 AxesLabel -> {"Restoring Force", "x[t]"}]

5970NonlinearSpring-revised.png
5970NonlinearSpring-revised.png360×226 2.4 KB

Note that in revising the notebook, WhenEvent is not needed. It is the physics, or if you like the math description of the physics, which I revised. Not the code.

Kind regards,
David

5

Now I understand it well. Thank you so much for the help David :slight_smile:
Regards, Khalid