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.sstatic.net/gN7PN.png

Improve this question
$\endgroup$

1 Answer 1

Reset to default
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.

Improve this answer
$\endgroup$
4
  • $\begingroup$ thanks you very much $\endgroup$ Commented Dec 26, 2020 at 11:42
  • $\begingroup$ You're welcome! Hope it helps $\endgroup$ Commented Dec 26, 2020 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$ Commented Dec 26, 2020 at 11:45
  • $\begingroup$ Evaluate a modified Table and use ListLinePlot3D[…] $\endgroup$ Commented Dec 26, 2020 at 11:49

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.