2
$\begingroup$

Questions regarding the use of nested For-loops have been posted a few times on StackExchange. However, I think that the solutions to those questions are not really applicable to what I am trying to do here. I am trying to use nested For-loops in NDSolve, which somehow does not work.

One of the methods included in the solution to other similar problems include things like putting If in Which. That does not seem to work for me.

My code:

ClearAll["Global`*"] z[t_] := y[t] - l Cos[φ[t]]; n[t_] := k Abs[z[t]]^(3/2) - c z'[t]; R = 0.0238; k = 141000000; c = 20; g = 9.81; μk = 0.1804; i = 8.92958*10^-6; m = 0.035; l = R Sqrt[2 (1 + Sin[(1/2) π - 2*0.41])]; a = NDSolve[ {m y''[t] == n[t] - m g, φ''[t] == (n[t]*l*Sin[φ[t]] + m x''[t]*l*Cos[φ[t]])/i, x''[t] == If[(x[t] - l*φ[t] == 0) // Evaluate, l (φ''[t] Cos[φ[t]] - φ'[t]^2 Sin[φ[t]]) // Evaluate, If[x[t] - l*φ[t] > 0, -μk*n[t]/m // Evaluate, -μk*n[t]/m // Evaluate]], y[0] == l*Cos[0.41], y'[0] == -2.22, φ[0] == 0.41, φ'[0] == -50, x[0] == 0, x'[0] == 0, WhenEvent[z[t] == 0 // Evaluate, "StopIntegration"; Print[y'[t]]; Print[φ'[t]]]}, {y'[t], y[t], x[t], x'[t], φ[t], φ'[t], z[t], n[t], φ''[t]}, {t, 0, 0.2}, Method -> {"EquationSimplification" -> "Residual", "DiscontinuityProcessing" -> False}, AccuracyGoal -> Automatic, WorkingPrecision -> MachinePrecision, MaxSteps -> 100000000, PrecisionGoal -> Automatic]

The piece of code containing the nested If is

x''[t] == If[(x[t] - l*φ[t] == 0) // Evaluate, l (φ''[t] Cos[φ[t]] - φ'[t]^2 Sin[φ[t]]) // Evaluate, If[x[t] - l*φ[t] > 0, -μk*n[t]/m // Evaluate, -μk*n[t]/m // Evaluate]],

It would be great if anyone can help.

Edit

I am trying to simulate a bounce of an object through this equation. I have the equation describing the motion along the y axis, x axis and rotation respectively. However, for the motion along the x-axis, it would be piecewise as during the collision, there may be both static friction and kinetic friction. I have to take both into account. Under kinetic friction, it will be either in the positive or negative direction depending on the relative motion between contact point and ground. That's why there are nested If expressions. The outer If is to differentiate between static and kinetic friction, while the inner If is to differentiate between positive and negative kinetic friction.

Improve this question
$\endgroup$
6
  • $\begingroup$ It would be helpful if before posting the code you formulate in words the problem you are going to solve with this code. If this is NDSolve, it should be a differential equation. Just formulate it. Describe, (either in words, or as a formula) why there should be conditional operators inside. $\endgroup$ Commented Dec 28, 2020 at 13:09
  • $\begingroup$ I have edited the post. Thanks $\endgroup$ Commented Dec 28, 2020 at 13:20
  • 1
    $\begingroup$ I don't think those calls to Evaluate are doing you any good. Read the Possible Issues section of the documentation article for Evaluate. $\endgroup$ Commented Dec 28, 2020 at 13:55
  • $\begingroup$ Have you seen Menu/Help/WolframDocumentation/NDSolve/Applications/Hybrid Differential Equations? $\endgroup$ Commented Dec 28, 2020 at 13:56
  • $\begingroup$ Yes I just took a look. Maybe WhenEvent is a good way, but I think it will come with a lot of limitation as well because lets say when the WhenEvent condition is satisfied, it will go from the first state to the second state. But What if the condition is satisfied again later on in the simulation, and i need mathematica to transit from the second state back to the first state. Can WhenEvent do that? @AlexeiBoulbitch $\endgroup$ Commented Dec 28, 2020 at 14:01

1 Answer 1

Reset to default
5
$\begingroup$

The problem is that z[t] and therefore n[t] depend on phi and y, but that dependence is internal and not explicit. E.g. the simplified code, something like which happens at the initial condition when NDSolve starts integrating, If[x[t] == 0, n[t], 2 n[t]] /. t -> 0 /. {x[0] -> 0, y[0] -> 0.0400369, \[CurlyPhi][0] -> 0.41} does not result in a number. The best solution is to rewrite z and n, e.g. z[y_, phi_] := y - l Cos[phi] and similarly for n[y_, phi_] := ...; then use n[y[t], phi[t]] in the equations. Then all the Evaluates are unnecessary, as @m_goldberg suggests in the comments. The quick fix is to put an Evaluate on the nested If:

NDSolve[{m y''[t] == n[t] - m g, φ''[ t] == (n[t]*l*Sin[φ[t]] + m x''[t]*l*Cos[φ[t]])/i, x''[t] == If[(x[t] - l*φ[t] == 0) // Evaluate, l (φ''[ t] Cos[φ[t]] - φ'[t]^2 Sin[φ[ t]]) // Evaluate, Evaluate@ If[x[t] - l*φ[t] > 0, -μk*n[t]/m // Evaluate, -μk*n[t]/m // Evaluate]], y[0] == l*Cos[0.41], y'[0] == -2.22, φ[0] == 0.41, φ'[0] == -50, x[0] == 0, x'[0] == 0, WhenEvent[z[t] == 0 // Evaluate, "StopIntegration"(*;Print[y'[t]]; Print[φ'[t]]*)]}, {y'[t], y[t], x[t], x'[t], φ[t], φ'[t], z[t], n[t], φ''[t]}, {t, 0, 0.2}, Method -> {"EquationSimplification" -> "Residual", "DiscontinuityProcessing" -> False}, AccuracyGoal -> Automatic, WorkingPrecision -> MachinePrecision, MaxSteps -> 100000000, PrecisionGoal -> Automatic]

Probably such a complicated problem has other issues. Certainly x[t] - l*φ[t] == 0 is probably never going to be True due to the discrete nature of numerical integration. It is given a special semantics in WhenEvent[x[t] - l*φ[t] == 0,...] but not in If[x[t] - l*φ[t] == 0, <one>, <two>]. The If variant is equivalent to simply <two> with a very high probability. So @bob_the_legend's comment concerning refactoring the Ifs with WhenEvent seems worth pursuing.

But to reiterate, I'd say the starting point is to follow good programming practice and make all parameters a function depends on be explicit arguments to the function.

Improve this answer
$\endgroup$
0

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.