1
$\begingroup$

It is impossible to solve numerically a nonlinear second-order equation, although the initial conditions are given.

ClearAll["Global`*"]

pars = {α = 0.15, h = 1, ω = 2 Pi 0.25, T = 1};

    extr = -(x[t] - 1)^2;
    
    func = extr + (D[
           x[t] - α Sin[ω t] + 
            D[x[t] - α Sin[ω t], t]/T, t] - 0)^2;
    
    sys1 = NDSolve[{x'[t] == 
         func α Sin[ω t] + D[α Sin[ω t], t], 
        x[0] == 2, x'[0] == -2}, x, {t, 0, 500}];

NDSolve::deqn: Equation or list of equations expected instead of False in the first argument {(x^′)[t]==0.235619 Cos[1.5708 t]+0.15 Sin[1.5708 t] (-Plus[<<2>>]^2+(Times[<<2>>]+Times[<<2>>]+(<<1>>^(<<1>>))[<<1>>]+(<<1>>^(<<1>>))[<<1>>])^2),x[0]==2,False}.
$\endgroup$
11
  • 1
    $\begingroup$ You need Clear[Derivative]. This isn't the only problem, though. $\endgroup$ – xzczd Dec 10 '20 at 9:15
  • 1
    $\begingroup$ One thing I can spot is, the ODE and the i.c. are inconsistent. Just try: {x'[t] == func \[Alpha] Sin[\[Omega] t] + D[\[Alpha] Sin[\[Omega] t], t], x[0] == 2, x'[0] == -2}/.t->0. If they should not, there's something wrong with your equation, please double check it. $\endgroup$ – xzczd Dec 10 '20 at 9:18
  • $\begingroup$ NDSolve::underdet: There are more dependent variables, {x[t]}, than equations, so the system is underdetermined. $\endgroup$ – dtn Dec 10 '20 at 9:20
  • 1
    $\begingroup$ I don't mean that you should feed the sample into NDSolve. Please execute it separately and observe the output, then you should understand what I mean. $\endgroup$ – xzczd Dec 10 '20 at 9:22
  • 1
    $\begingroup$ Once again, please observe the output. $\endgroup$ – xzczd Dec 10 '20 at 9:24
1
$\begingroup$

Here, we take care of the singularity at t == 0 and let "StiffnessSwitching" + a high working precision slide past the rest of them, which occur whenever t is equal to an even integer.

Set up:

Block[{α = 15/100, h = 1, ω = 2 Pi /4, T = 1},
 extr = -(x[t] - 1)^2;
 func = extr + (D[
       x[t] - α Sin[ω t] + 
        D[x[t] - α Sin[ω t], t]/T, t] - 0)^2;
 
 ode = {x'[t] == 
    func α Sin[ω t] + D[α Sin[ω t], t]};
 ics = {x[0] == 2, x'[0] == -2};
 ]

Form a solution of the form $x(t) = u(t) + \Delta x(t)$, where $u''(0)$ is finite and $\Delta x''(0)$ is infinite. Then derive the ODE for $u$, integrate, and form the solution for $x$:

xpp = Solve[ode, x''[t]];
xpp /. Thread[{x[t], 
     x'[t]} -> ({x[t], x'[t]} /. t -> 0 /. Solve@ics // First)];
% // Simplify;
xpp0 = x''[t] /. %;

dx = (Series[#, {t, 0, 0}, Assumptions -> t > 0] & /@ 
     xpp0) /.
   {HoldPattern[
      SeriesData[t_, t0_, coeff_, min_?Negative, max_, denom_]] :>
     Total@Take[coeff*(t - t0)^(Range[min, max - 1]/denom), -min],
    _SeriesData :> 0};
dx = Integrate[dx, t, t]

{-(4/3) I Sqrt[1 + 80/(3 π)] t^(3/2), 
 4/3 I Sqrt[1 + 80/(3 π)] t^(3/2)}

uodes = MapThread[
    Function[{pp, dx}, 
     x''[t] == pp /. {x -> (u[#] + dx &)}], {x''[t] /. 
      xpp, (dx /. t -> #)}
    ] // Simplify;
upp = First@Solve[#, u''[t]] & /@ uodes;

uppic = MapThread[
  Function[{pp, dx},
   Limit[u''[t] /. pp /. 
     Thread[{u[t], 
        u'[t]} -> ({u[0], u'[0]} /. 
         First@Solve[ics /. x -> (u[#] + dx &)])],
    t -> 0,
    Direction -> "FromAbove"]
   ],
  {upp, dx /. t -> #}]
(* only one branch satisfies IC *)

(* {DirectedInfinity[I], 2 + 40/(3*Pi) + (3*Pi)/40}  *)

obj[t_?NumericQ, u_?NumericQ, up_?NumericQ] =
  Piecewise[
    {{u''[t] /. Last@upp, t > 0}},
    Last@uppic] /. {u[t] -> u, u'[t] -> up};
uics = Equal @@@ 
   First@Solve[ics /. x -> (Evaluate[u[#] + Last@dx /. t -> #] &)];

usol = NDSolve[{u''[t] == obj[t, u[t], u'[t]], uics},
    u, {t, 0, 100.1},
    InterpolationOrder -> All,
    Method -> {"StiffnessSwitching", 
      Method -> {"ExplicitRungeKutta", Automatic}},
    MaxSteps -> 100000,
    "ExtrapolationHandler" -> {Indeterminate &, 
      "WarningMessage" -> False}
    , WorkingPrecision -> 30, PrecisionGoal -> 8
    ]; // AbsoluteTiming

(*  {44.0669, Null}  *)

xsol = {x -> (Evaluate[u[#] + Last@dx /. t -> #] &)} /. usol;

Check initial conditions:

ics /. xsol

(*  {{True, True}}  *)

Visualize:

ReImPlot[x[t] /. xsol // Evaluate, {t, 0, 100}]

enter image description here

Residual error on the original ODE:

Plot[ode /. Equal -> Subtract /. xsol // RealExponent // Evaluate,
 {t, 0, 100}]

enter image description here

$\endgroup$
1
$\begingroup$

To long for a comment: The initial conditions aren't consistent!

The ode at t==0 evaluates to

x'[t] == func \[Alpha] Sin[\[Omega] t] + D[\[Alpha] Sin[\[Omega] t], t] /. t -> 0  
(*Derivative[1][x][0] == 0.235619*)  

which differs from your second initial condition x'[0]==2!

Be aware, that your ode switches between first and second order. That seems to be hard to handle by NDSolve

Some of the error messages NDSolve might be bypassed if you use Method -> {"EquationSimplification" -> "Residual"}:

X = NDSolveValue[{x'[t] ==func \[Alpha] Sin[\[Omega] t] + D[\[Alpha] Sin[\[Omega] t], t],x[0] == 2, x'[0] == -2}, x, {t, 0, 1}, 
Method -> {"EquationSimplification" -> "Residual"}]   
$\endgroup$
0

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.