2
$\begingroup$

I am trying to simulate something using differential equations. However, the equations are piecewise. Here is my code:

z[t_] = y[t] - l Cos[φ[t]];
n[t_] = k Abs[z[t]] - k2 z'[t];
k = 50000000;
k2 = 37;
g = 9.81;
μ = 0.3;
i = 10^-5;
l = 0.01;
m = 0.35;
a = NDSolve[
  If[μ >= (m x''[t])/n[t], (*Condition*)
(*First set of differential equation*)  
{m y''[t] == n[t] - m g,
    φ''[t] == (
     n[t] l Sin[φ[t]] + m x''[t] Cos[φ[t]])/i,
    x''[t] == 
     l (φ''[t] Cos[φ[t]] - φ'[
          t]^2 Sin[φ[t]])},
(*second set of differential equation*) 
   {φ''[t] == (
     n[t] l Sin[φ[t]] + m x''[t] Cos[φ[t]])/i,
    m y''[t] == n[t] - m g,
    m x''[t] == μ n[t]}],
  
  {y[0] == 0.01, 
   y'[0] == -3, φ[0] == 0.01, φ'[0] == 0, 
   x[0] == 0, x'[0] == 0},
  {y[t], φ[t], x[t]}, {t, 0, 0.0009}]

I am trying to do a numerical simulation:

  1. When the If condition is satisfied, I would like Mathematica to solve the first set of differential equations.
  2. When the If conditions is no longer satisfied, I would like Mathematica to continue the simulation with the second set of differential equations.
  3. When the If condition is satisfied again, it will continue the simulation with the first set of differential equations again. The code will stop running when n[t] becomes 0 (which happens to be the initial conditions as well) .
  4. Last but not least, I want to plot the graphs of the functions.

I have checked through the two sets of differential equations individually using NDSolve and they are all solvable.

I am relatively new to Mathematica and still learning. It would be great if anyone can help. Thanks.

$\endgroup$
2
  • 2
    $\begingroup$ WhenEvent maybe one approach. $\endgroup$
    – cvgmt
    Dec 21, 2020 at 4:55
  • $\begingroup$ Thanks, I will read up more on that $\endgroup$ Dec 21, 2020 at 4:56

1 Answer 1

2
$\begingroup$

Thanks your question. I am also the first time use WhenEvent and If in ODE,so I am not sure whether it is right or not.

At the beginning,I consider WhenEvent,however I cann't find the way. I only get the time t when the event μ*n[t] >= m x''[t] occure. I using the code WhenEvent[μ*n[t] >= m x''[t] // Evaluate, Print[t]]

z[t_] = y[t] - l Cos[φ[t]];
n[t_] = k Abs[z[t]] - k2 z'[t];
k = 50000000;
k2 = 37;
g = 9.81;
μ = 0.3;
i = 10^-5;
l = 0.01;
m = 0.35;
a = NDSolveValue[{m y''[t] == 
    n[t] - m g, φ''[
     t] == (n[t] l Sin[φ[t]] + 
       m x''[t] Cos[φ[t]])/i, 
   x''[t] == l (φ''[t] Cos[φ[t]] - φ'[t]^2 Sin[φ[t]]), 
y[0] == 0.01, y'[0] == -3, φ[0] == 0.01, φ'[0] == 0, 
   x[0] == 0, x'[0] == 0, 
   WhenEvent[μ*n[t] >= m x''[t] // Evaluate, Print[t]]}, {y[
    t], φ[t], x[t]}, {t, 0, 0.0009}]

(* 0.000263796 *)

Then I back to If and change the equation m*x''[t]==μ n[t] to x''[t]==μ n[t]/m,as below. The main settings is

  x''[t] == 
    If[(μ*n[t] >= m x''[t]) // Evaluate, (μ n[t]/m) //Evaluate, 
     l (φ''[t] Cos[φ[t]] - φ'[t]^2 Sin[φ[t]])]
z[t_] := y[t] - l Cos[φ[t]];
n[t_] := k Abs[z[t]] - k2 z'[t];
k = 50000000;
k2 = 37;
g = 9.81;
μ = 0.3;
i = 10^-5;
l = 0.01;
m = 0.35;
a = NDSolveValue[{m y''[t] == 
    n[t] - m g, φ''[
     t] == (n[t] l Sin[φ[t]] + 
       m x''[t] Cos[φ[t]])/i, 
   x''[t] == 
    If[(μ*n[t] >= m x''[t]) // Evaluate, (μ n[t]/m) // 
      Evaluate, 
     l (φ''[t] Cos[φ[t]] - φ'[t]^2 Sin[φ[t]])],
 y[0] == 0.01, y'[0] == -3, φ[0] == 0.01, φ'[0] == 0, 
   x[0] == 0, x'[0] == 0}, {y[t], φ[t], x[t]}, {t, 0, 
   0.0009}, 
  Method -> {"EquationSimplification" -> "Residual", 
    "DiscontinuityProcessing" -> False}]
Plot[%, {t, 0, 0.0009}]
$\endgroup$
12
  • $\begingroup$ Hi,thanks for your response. However, when I run it on my mathematica, there are some error messages. Also may i know why is it 3Dplot? I would like to have 3 different graphs of the 3 functions. thank you :) $\endgroup$ Dec 21, 2020 at 12:03
  • $\begingroup$ Also, I would like to know if it is possible to stop the graph when n[t] becomes 0? Thanks $\endgroup$ Dec 21, 2020 at 12:10
  • $\begingroup$ @bobthelegend I condiser n[t] becomes 0 later. I think at that time we can use WhenEvent $\endgroup$
    – cvgmt
    Dec 21, 2020 at 12:14
  • $\begingroup$ Okay i will try to do that thanks $\endgroup$ Dec 21, 2020 at 12:18
  • 1
    $\begingroup$ @bobthelegend Maybe we can add WhenEvent[n[t] == 0 // Evaluate, "StopIntegration"] to stop $\endgroup$
    – cvgmt
    Dec 21, 2020 at 15:44

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.