NDSolve::ndnum for initial value problem of 2nd order ODE

Question

I have this code, which I think should work fine

Param = {Mo -> 112, mu1 -> 275.52,k -> 4.144,m2 -> 14,m3 -> 24, m4 -> 12, 
         R3 -> 0.16, r3 -> 0.12, rho3 -> 0.14, r1 -> 0.06};
eq = (m3*r1^2 + m3 *(r1*(r3+R3)/(2*r3+R3)) ^2 + m3*rho3*(r1/(2*r3+R3))^2 + 
      m4*(r1*(r3+R3)/(2*r3+R3))^2)*(phi1''[t]*Sin[phi1[t]]^2 - 
      0.5*phi1'[t]*Sin[phi1[t]]*Cos[phi1[t]]);
Q = Mo - k*phi1'[t] + mu1*phi1'[t]*4*(r1*(r3+R3)/(2*r3+R3));
sol = NDSolve [{eq==Q, phi1[0] == 0, phi1'[0] == 23.375} /.Param,
               {phi1},{t, 0,0.2688}] 

However, it gives me warning

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`.

There are a lot of similar questions, but none of them were helpful since I don't understand what that error means and what causes it. I believe it has something to do with phi1[t] being inside a sine function.

Answer 1

Let see what happens to the ode at $t=0$

In your case

ClearAll["Global`*"]
Param = {Mo -> 112, mu1 -> 275.52, k -> 4.144, m2 -> 14, m3 -> 24, 
  m4 -> 12, R3 -> 0.16, r3 -> 0.12, rho3 -> 0.14, r1 -> 0.06}
eq = (m3*r1^2 + m3*(r1*(r3 + R3)/(2*r3 + R3))^2 + 
    m3*rho3*(r1/(2*r3 + R3))^2 + 
    m4*(r1*(r3 + R3)/(2*r3 + R3))^2)*(phi1''[t]*Sin[phi1[t]]^2 - 
    0.5*phi1'[t]*Sin[phi1[t]]*Cos[phi1[t]])
Q = Mo - k*phi1'[t] + mu1*phi1'[t]*4*(r1*(r3 + R3)/(2*r3 + R3));
ode = eq == Q 
ode = ode /. Param // Rationalize
ode = First@ode - Last@ode == 0

You can see that at $t=0$ the above reduces to

Because phi1[t]=0 at $t=0$, and hence $\sin(0)=0$ so all of the last term goes away.

Now lets apply the other initial condition which is phi1'[0]== 187/8 and see if it satisfies the above at $t=0$

ode /. phi1[t] -> 0

% /. phi1'[t] -> Rationalize[23.375]

So at $t=0$, the initial conditions are not consistent with the ode.

So I think the problem is that your initial conditions are not consistent. Normally NDSolve will say that as another warning message ibcinc , but may be in this case it did not?

Answer 2

The problem is strongly related to e.g.

Numerical solution of Bessel-like equation using NDSolve

NDSolve divide-by-zero trouble

Trouble in the second order ordinary differential equations with second oder coupled iterm

But it may be hard for beginners to notice the underlying relationship, so let me elaborate a bit. As discussed in this post, NDSolve by default transforms the ODE to a first order system in standard form. This standard form can be obtained with to1storder:

internalsys = 
 to1storder[eq == Q /. Param, phi1, t, "standardize" -> True, "form" -> Subscript]

As shown above, the standard form involves a Csc[Subscript[phi1, 0][t]] term. (Subscript[phi1, 0] stands for phi1 in original equation. ) So the initial condition phi1[0] == 0 will lead to non-numerical ComplexInfinity:

internalsys[[1]] /. t -> 0 /. Subscript[phi1, 0][0] -> 0
(* {Subscript[phi1, 1]'[0] == ComplexInfinity} *)

This is the reason why ndnum warning pops up.

The ndnum warning can be circumvented by, as discussed in linked posts above, SolveDelayed -> True:

NDSolve[{eq == Q, phi1[0] == 0, phi1'[0] == 23.375} /. Param, {phi1}, {t, 0, 0.2688},
  SolveDelayed -> True]

But this leads to the warning

NDSolve::ivres: NDSolve has computed initial values that give a zero residual for the differential-algebraic system, but some components are different from those specified. If you need them to be satisfied, giving initial conditions for all dependent variables and their derivatives is recommended.

The reason is already elaborated in Nasser's answer, so I'd like not to repeat. You should double check the underlying model of your equation.