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 '20 at 4:55
  • $\begingroup$ Thanks, I will read up more on that $\endgroup$ – bob the legend Dec 21 '20 at 4:56
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$ – bob the legend Dec 21 '20 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$ – bob the legend Dec 21 '20 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 '20 at 12:14
  • $\begingroup$ Okay i will try to do that thanks $\endgroup$ – bob the legend Dec 21 '20 at 12:18
  • 1
    $\begingroup$ @bobthelegend Maybe we can add WhenEvent[n[t] == 0 // Evaluate, "StopIntegration"] to stop $\endgroup$ – cvgmt Dec 21 '20 at 15:44

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.