4
$\begingroup$

I am looking for advice on the use of NDsolve in the case where there are two coupled ODEs with some derivatives not being written explicitly. The ones that aren't written explicitly occur within a square root function. In that form I cannot seem to get NDSolve to run. I would be grateful for any assistance or advice.

To give some background, I modelled a novel type of catapult in Mathematica using Lagrangian dynamics.The model did not account for air drag losses. The type uses a spinning high inertia rotor with a projectile attached by a sling to the end of the rotor. The slinging out of the projectile draws energy from the spinning rotor similarly to a trebuchet.

I then built an actual machine and compared the performance with the model. In doing so it is clear that the actual trajectory lags behind the predicted trajectory due to air drag not being modelled. Note that I am not referring to the trajectory after the projectile is launched; I refer to the trajectory from time t=o until the program stops at the time of the launch.

I am trying model air drag now and can successfully incorporate linear drag into the model. However for the quantities under consideration ( 2 kg sphere of diameter 160 mm launched at 100 meters/second) quadratic drag will dominate.

It is the modelling of quadratic drag that requires a radical containing derivatives. That is where I am having difficulty with NDSolve. The term in question is the norm of the velocity vector and would appear in both equations as Sqrt[(a1*Theta'[t])^2+(s*Phi'[t])^2+Cos[Theta[t]-Phi[t]]*Theta'[t]*Phi'[t]]. I believe it is the derivatives occurring inside the radical sign that are causing the issues but I am unable to simplify the equations or put them in a form that NDSolve can process. The code modelling the linear drag is shown below and the equations I am having difficulty with are shown below that. Thank you for any advice that may be offered.

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Manipulate[
 a = 2.035;                                     "Arm Length";
 s = 2.261;                                     "Sling Length";
 i = 200;                                         "Arm Moment of \
Inertia";
 m = 2;                                             "Projectile Mass";
 sia = 3.5925 + 
   aia ;                 "Sling Initial Angle and Arm Initial Angle";

 eqns = {θ''[
     t] == (-m*a*s*Cos[θ[t] - ϕ[t]]/(i + m*a^2)) ϕ''[
       t] + ( -m*a*s*
        Sin[θ[t] - ϕ[t]]/(i + m*a^2)) ϕ'[t]^2 - 
     k (a*s*Cos[θ[t] - ϕ[t]]/(i + m*a^2)*ϕ'[t] + 
        a^2/(i + m*a^2) θ'[t]),
              ϕ''[
     t] == (-a*Cos[θ[t] - ϕ[t]]/s ) θ''[
       t] + (a*Sin[θ[t] - ϕ[t]]/s)*θ'[
        t]^2                                                          \
              - 
     k (a*Cos[θ[t] - ϕ[t]]/(m*s)*θ'[
          t] + (1/m) ϕ'[
          t])                                        };

 sol = NDSolve[{eqns,
    θ[0] == aia, ϕ[0] == sia, θ'[0] == 
     aiv, ϕ'[0] == aiv},
     {θ, ϕ}, {t, 5},
     Method -> {"EventLocator", 
     "Event" -> (ϕ[t] - θ[t]) - (Pi + β), 
     "EventAction" :> Throw[tend = t, "StopIntegration"]}];

 x[t_] := a*Cos[θ[t]] + 
   s*Cos[ϕ[t]];                       "x Coordinate";
 y[t_] := 
  a*Sin[θ[t]] + 
   s*Sin[ϕ[t]];                         "y Coordinate";
 z[t_] := {x[t], 
   y[t]};                                                        \
"Position";
 vel[t_] := 
  Sqrt[x'[t]^2 + 
    y'[t]^2];                           "Projectile Speed";
 accel[t_] := 
  Sqrt[x''[t]^2 + y''[t]^2]/g;              "Projectile Accelleration";

 frames = Table[Circle[z[ff + (i - 1)*.03333], .05], {i, f}] /. sol;
 box = Line[{{-4.3, -4.3}, {-4.3, 4.3}, {4.3, 
     4.3}, {4.3, -4.3}, {-4.3, -4.3}}];
 trajectory = 
  ParametricPlot[z[t] /. sol, {t, 0, tend}, PlotStyle -> Red][[1]];
 sling = Line[{{a1*Cos[θ[tend]], a*Sin[θ[tend]]}, 
     z[tend]}] /. sol;
 arm = Evaluate[
    Line[{{0, 0}, {a*Cos[θ[tend]], 
       a*Sin[θ[tend]]}}]] /. sol;

 Grid[{{"Launch Elev. (Deg)", 
     Round[-(180/Pi)*Evaluate[ArcTan[y'[tend]/x'[tend]]] /. sol, .1][[
      1]]},
    {"Acceleration (g's)", Round[accel[tend] /. sol, 1][[1]]},
    {"Velocity (m/sec)" , Round[vel[tend] /. sol, .1][[1]]},
    {"Sling Tension (Kg)", 
     Round[Round[m*accel[tend] /. sol, 1][[1]], .1]},
    {"θ'Final (rad/sec)", 
     Round[Evaluate[θ'[tend]] /. sol, .1][[1]]},
    {"ϕ'Final(rad/sec) ", 
     Round[Evaluate[ϕ'[tend]] /. sol, .1][[1]]}, 
    {"Initial Energy (kJ)", 
     Round[Evaluate[.0005*i*(θ'[0])^2] /. sol, .1][[1]]},
    {"Elapsed Time", Round[tend, .001]}},
   Frame -> All, ItemSize -> 11]



  Graphics[{box, frames, trajectory, {Green, sling}, {Thick, arm}}, 
   PlotRange -> {{-5, 5}, {-5, 5}}, ImageSize -> 450],
 {{k, 0, " Linear Drag Coefficient"}, 0, 15, .01, 
  Appearance -> "Labeled"},
 {{ff, 0, "First Frame"}, 0, 0.0333, .0001, Appearance -> "Labeled"},
 {{f, 1, "Number of Frames"}, 1, 10, 1, Appearance -> "Labeled"},
 Delimiter,
 {{aia, 4.86, "Initial Position"}, 0, 2*Pi, .01, 
  Appearance -> "Labeled"},
 {{aiv, 17.78, "Initial Angular Velocity "}, .1, 40, .01, 
  Appearance -> "Labeled"},
 {{β, 2.7, "Release Angle"}, 0, 2 Pi, .05, 
  Appearance -> "Labeled"},
 TrackedSymbols -> Manipulate, ContinuousAction -> True]

The equations that I am having difficulty with are :

eqns = {θ''[
     t] == (-m*a*s*Cos[θ[t] - ϕ[t]]/(i + m*a^2)) ϕ''[
       t] + ( -m*a*s*
        Sin[θ[t] - ϕ[t]]/(i + m*a^2)) ϕ'[t]^2 - 
     k (Sqrt[a^2*θ'[t]^2 + s^2*ϕ'[t]^2 + 
         Cos[θ[t] - ϕ[t]]*θ'[t]*ϕ'[t]]) (a*s*
         Cos[θ[t] - ϕ[t]]/(i + m*a^2)*ϕ'[t] + 
        a^2/(i + m*a^2) θ'[t]),
              ϕ''[
     t] == (-a*Cos[θ[t] - ϕ[t]]/s ) θ''[
       t] + (a*Sin[θ[t] - ϕ[t]]/s)*θ'[
        t]^2                                                          \
         - 
     k (Sqrt[a^2*θ'[t]^2 + s^2*ϕ'[t]^2 + 
         Cos[θ[t] - ϕ[t]]*θ'[t]*ϕ'[t]]) (a*
         Cos[θ[t] - ϕ[t]]/(m*s)*θ'[
          t] + (1/m) ϕ'[
          t])                                        };
$\endgroup$
  • $\begingroup$ Could you reduce this to a minimal working example? $\endgroup$ – Sascha Sep 9 '16 at 14:56
  • $\begingroup$ Two typos in your Manipulate. You are multiplying the Grid by the Graphics. Enclose them in a Column. You have a1 in the equation for sling. It should be a. $\endgroup$ – Jack LaVigne Sep 9 '16 at 15:16
  • $\begingroup$ If your drag is quadratic, why do you use radical? $\endgroup$ – Vsevolod A. Mar 15 '18 at 14:07
3
$\begingroup$

I am Mathematica 11.0 and windows 7.

I copied your Manipulate and fixed the two typos mentioned in my comment.

I copied your equation with drag and pasted it into the Manipulate.

I did not load the package Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];

It appeared to work fine in the sense that NDSolve spit out an answer.

Manipulate[
 a = 2.035;(* Arm Length *)
 s = 2.261;(* Sling Length *)
 i = 200;(* Arm Moment of Inertia *)
 m = 2;(* Projectile Mass *)
 sia = 3.5925 + aia;(* Sling Initial Angle and Arm Initial Angle *)

 eqns = {θ''[
     t] == (-m*a*s*Cos[θ[t] - ϕ[t]]/(i + m*a^2)) ϕ''[
       t] + (-m*a*s*Sin[θ[t] - ϕ[t]]/(i + m*a^2)) ϕ'[
        t]^2 - k (Sqrt[
        a^2*θ'[t]^2 + s^2*ϕ'[t]^2 + 
         Cos[θ[t] - ϕ[t]]*θ'[t]*ϕ'[t]]) (a*s*
         Cos[θ[t] - ϕ[t]]/(i + m*a^2)*ϕ'[t] + 
        a^2/(i + m*a^2) θ'[t]), ϕ''[
     t] == (-a*Cos[θ[t] - ϕ[t]]/s) θ''[
       t] + (a*Sin[θ[t] - ϕ[t]]/s)*θ'[t]^2 - 
     k (Sqrt[a^2*θ'[t]^2 + s^2*ϕ'[t]^2 + 
         Cos[θ[t] - ϕ[t]]*θ'[t]*ϕ'[t]]) (a*
         Cos[θ[t] - ϕ[t]]/(m*s)*θ'[
          t] + (1/m) ϕ'[t])};

 sol = NDSolve[{eqns, θ[0] == aia, ϕ[0] == 
     sia, θ'[0] == aiv, ϕ'[0] == 
     aiv}, {θ, ϕ}, {t, 5}, 
   Method -> {"EventLocator", 
     "Event" -> (ϕ[t] - θ[t]) - (Pi + β), 
     "EventAction" :> Throw[tend = t, "StopIntegration"]}];

 x[t_] := a*Cos[θ[t]] + s*Cos[ϕ[t]]; "x Coordinate";
 y[t_] := a*Sin[θ[t]] + s*Sin[ϕ[t]]; "y Coordinate";
 z[t_] := {x[t], y[t]}; "Position";
 vel[t_] := Sqrt[x'[t]^2 + y'[t]^2]; "Projectile Speed";
 accel[t_] := Sqrt[x''[t]^2 + y''[t]^2]/g; "Projectile Accelleration";
 frames = Table[Circle[z[ff + (i - 1)*.03333], .05], {i, f}] /. sol;
 box = Line[{{-4.3, -4.3}, {-4.3, 4.3}, {4.3, 
     4.3}, {4.3, -4.3}, {-4.3, -4.3}}];
 trajectory = 
  ParametricPlot[z[t] /. sol, {t, 0, tend}, PlotStyle -> Red][[1]];
 sling = Line[{{a*Cos[θ[tend]], a*Sin[θ[tend]]}, 
     z[tend]}] /. sol;
 arm = Evaluate[
    Line[{{0, 0}, {a*Cos[θ[tend]], 
       a*Sin[θ[tend]]}}]] /. sol;

 Column[{
   Grid[
    {
     {"Launch Elev. (Deg)",
      Round[-(180/Pi)*Evaluate[ArcTan[y'[tend]/x'[tend]]] /. 
         sol, .1][[1]]},
     {"Acceleration (g's)", Round[accel[tend] /. sol, 1][[1]]},
     {"Velocity (m/sec)", Round[vel[tend] /. sol, .1][[1]]},
     {"Sling Tension (Kg)", 
      Round[Round[m*accel[tend] /. sol, 1][[1]], .1]},
     {"θ'Final (rad/sec)", 
      Round[Evaluate[θ'[tend]] /. sol, .1][[1]]},
     {"ϕ'Final(rad/sec) ", 
      Round[Evaluate[ϕ'[tend]] /. sol, .1][[1]]},
     {"Initial Energy (kJ)", 
      Round[Evaluate[.0005*i*(θ'[0])^2] /. sol, .1][[1]]},
     {"Elapsed Time", Round[tend, .001]}
     },
    Frame -> All,
    ItemSize -> 11
    ],

   Graphics[
    {box,
     frames,
     trajectory,
     {Green, sling},
     {Thick, arm}},
    PlotRange -> {{-5, 5}, {-5, 5}},
    ImageSize -> 450]
   }],

 {{k, 1, " Linear Drag Coefficient"}, 0, 15, .01, 
  Appearance -> "Labeled"},
 {{ff, 0, "First Frame"}, 0, 0.0333, .0001, Appearance -> "Labeled"},
 {{f, 1, "Number of Frames"}, 1, 10, 1, Appearance -> "Labeled"},
 Delimiter,
 {{aia, 4.86, "Initial Position"}, 0, 2*Pi, .01, 
  Appearance -> "Labeled"},
 {{aiv, 17.78, "Initial Angular Velocity "}, .1, 40, .01, 
  Appearance -> "Labeled"},
 {{β, 2.7, "Release Angle"}, 0, 2 Pi, .05, 
  Appearance -> "Labeled"},

 TrackedSymbols -> Manipulate,
 ContinuousAction -> True
 ]

Mathematica graphics

$\endgroup$
  • $\begingroup$ That works well now and is quite a bit tidier. Thank you very much, Jack, for your assistance and improvements. $\endgroup$ – JCF Sep 10 '16 at 2:04

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.