2
$\begingroup$

I want to plot a parametric 3D plot from the following ParametricNDSolve code:

ClearAll["Global`*"]
z[t_] := y[t] - l Cos[φ[t]];
n[t_] := k Abs[z[t]]^(3/2) - c z'[t];
R = 0.01395;
k = 1500000000;
c = 0;
g = 9.81;
μs = 0.3;
μk = 0.301;
i = 8.784*10^-6;
m = 0.035;
l = R Sqrt[2 (1 + Sin[1/2 π - 2*f])];
a = ParametricNDSolve[
  {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, μk *n[t]/m // Evaluate],
   y[0] == l*Cos[f], 
   y'[0] == -2.9, φ[0] == f, φ'[0] == h, 
   x[0] == 0, x'[0] == 0, 
   WhenEvent[z[t] == 0 // Evaluate ,end=t; "StopIntegration";]}, {y'[
    t]}, {t, 0, 0.01}, {f, h}, 
  Method -> {"EquationSimplification" -> "Residual", 
    "DiscontinuityProcessing" -> False}, AccuracyGoal -> Automatic, 
  WorkingPrecision -> MachinePrecision, MaxSteps -> 100000000, 
  PrecisionGoal -> Automatic]

The code works. However, when I tried to plot y'[end] against f and h using the following code: ParametricPlot3D[y'[end], {f, 0, 0.7}, {h, 0, 10}]

It outputs an empty graph without any error message. [![enter image description here][1]][1] I think the problem is because of the end=t in my WhenEvent. However, I have no idea how to fix it.

Please send help. Thanks a lot. [1]: https://i.stack.imgur.com/gN7PN.png

$\endgroup$
4
$\begingroup$

If you only need y'[end] as a result modify the ParametricNDSolve to

ysend = ParametricNDSolveValue[{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, μk*n[t]/m// Evaluate], 
y[0] == l*Cos[f],y'[0] == -2.9, φ[0] == f, φ'[0] == h,x[0] == 0, x'[0] == 0,WhenEvent[z[t] == 0 // Evaluate, end = t; "StopIntegration";]}, 
y'[end], {t, 0, 0.01}, {f, h},Method -> {"EquationSimplification" -> "Residual","DiscontinuityProcessing" -> False}, AccuracyGoal -> Automatic,WorkingPrecision -> MachinePrecision, MaxSteps -> 100000000,PrecisionGoal -> Automatic]

Now ysend[f, h] corresponds to y'[end] .

The evaluation is very time consuming, that's why I only give a table of selected results

Table[{f, h, ysend[f, h] }, {f, Subdivide[0, 0.7, 3]}, {h,Subdivide[0, 10, 3]}] // Quiet
(*{{{0., 0, 2.89997}, {0., 10/3, 2.87237}, {0., 20/3, 2.87005}, 
{0., 10,2.86764}}, {{0.233333, 0, 1.35025}, {0.233333, 10/3,1.31775},
{0.233333, 20/3, 1.28522}, {0.233333, 10,1.25266}}, {{0.466667, 0, 0.307301},
{0.466667, 10/3,0.264995}, {0.466667, 20/3, 0.222678}, {0.466667, 10,0.180351}}, 
{{0.7, 0, -0.0400896}, {0.7, 10/3, -0.0858449}, {0.7, 20/3, -0.131606}, {0.7, 10, -0.177372}}}*)

Plot3D[ ysend[f, h] , {f, 0, 0.7}, {h, 0, 10}] // Quiet should work but needs endless simulation time.

$\endgroup$
4
  • $\begingroup$ thanks you very much $\endgroup$ Dec 26 '20 at 11:42
  • $\begingroup$ You're welcome! Hope it helps $\endgroup$ Dec 26 '20 at 11:42
  • $\begingroup$ Is there anyway to make the simulation shorter. I dont need to get a very smooth 3D plot. a rough one showing the general trend would do. :) $\endgroup$ Dec 26 '20 at 11:45
  • $\begingroup$ Evaluate a modified Table and use ListLinePlot3D[…] $\endgroup$ Dec 26 '20 at 11:49

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.