Why does NDSolve[] want to treat my equations as Delay Differential Equations?

I use NDSolve command (method->equation simplification -> residual) to solve different sets of Ordinary differential equations (ODEs).

It had worked alright till I added rather a simple component Vspr. From then NDSolve[] cant solve such ODEs - I get an error: *NDSolve: Initial history needs to be specified for all variables for delay-differential equations.

I'm in fact sure that my equations didnt turn to be delay-differential equations. It must be an internal fail of NDSolve[] solver.

I've read two posts on this forum which concerns such a problem - NDSolve error: Cannot solve to find an explicit formula for the derivatives. Consider using the option Method->{"EquationSimplification"->"Residual"} and Solving this equation with NDSolve. The answers to them was to get rid of integrating inside ODEs. I tried somewhat to apply it to my problem nevertheless I dont have explicit integration in my ODEs.

Needs["Quaternions"]
angleVector[modul_, angleOYminus_] := {modul Cos[angleOYminus - Pi/2],
modul Sin[angleOYminus - Pi/2], 0}
modul[vector_] :=
Sqrt[(vector[[1]])^2 + (vector[[2]])^2 + (vector[[3]])^2]
(*kinematics:*)
tauP93third[l1_, rbp_, xbp_, ybp_] := \[Pi] -
ArcCsc[Sqrt[xbp^2 + ybp^2]/(l1 - rbp)] +
ArcSin[ybp/Sqrt[xbp^2 + ybp^2]]
łP10P9start = ls0 + ll2;
yP1st21sixth[t_] := -ll2 - ls0 +
rbp (1/2 (\[Pi] + 2 ArcSin[ybp/Sqrt[xbp^2 + ybp^2]] +
2 ArcSin[
rbp/(\[Sqrt]((-ybp -
l1 Cos[tauP10st1st21plus[t]])^2 + (-xbp +
l1 Sin[tauP10st1st21plus[t]])^2))] +
2 ArcSin[(l1 Sin[
ArcSin[ybp/Sqrt[xbp^2 + ybp^2]] +
1/2 (\[Pi] - 2 tauP10st1st21plus[t])])/(\[Sqrt]((-ybp -
l1 Cos[tauP10st1st21plus[t]])^2 + (-xbp +
l1 Sin[tauP10st1st21plus[t]])^2))]) +
1/2 (-\[Pi] - 2 ArcSin[ybp/Sqrt[xbp^2 + ybp^2]] -
2 ArcSin[
rbp/(\[Sqrt]((-ybp -
l1 Cos[tauP10st1st21plus[t2]])^2 + (-xbp +
l1 Sin[tauP10st1st21plus[t2]])^2))] -
2 ArcSin[(l1 Sin[
ArcSin[ybp/Sqrt[xbp^2 + ybp^2]] +
1/2 (\[Pi] -
2 tauP10st1st21plus[t2])])/(\[Sqrt]((-ybp -
l1 Cos[tauP10st1st21plus[t2]])^2 + (-xbp +
l1 Sin[tauP10st1st21plus[
t2]])^2))])) - \[Sqrt]((-ybp +
l1 Cos[ArcSin[(l1 - rbp)/Sqrt[xbp^2 + ybp^2]] -
ArcSin[ybp/Sqrt[xbp^2 + ybp^2]]] -
rbp Cos[ArcSin[(l1 - rbp)/Sqrt[xbp^2 + ybp^2]] -
ArcSin[ybp/Sqrt[xbp^2 + ybp^2]]])^2 + (-xbp +
l1 Sin[ArcSin[(l1 - rbp)/Sqrt[xbp^2 + ybp^2]] -
ArcSin[ybp/Sqrt[xbp^2 + ybp^2]]] -
rbp Sin[ArcSin[(l1 - rbp)/Sqrt[xbp^2 + ybp^2]] -
ArcSin[ybp/Sqrt[xbp^2 + ybp^2]]])^2) + \[Sqrt](-rbp^2 +
xbp^2 + ybp^2 + 2 l1 ybp Cos[tauP10st1st21plus[t]] +
l1^2 Cos[tauP10st1st21plus[t]]^2 -
2 l1 xbp Sin[tauP10st1st21plus[t]] +
l1^2 Sin[tauP10st1st21plus[t]]^2)
P1[t_] := {xbp + rbp, ybp - ybpP1tstart + yP1st21sixth[t], 0}
tauw[t_] := {0, 0, tauP10st1st21plus[t]}

versorz = UnitVector[3, 3];
skalartau = versorz.tauw[t];

l1w = angleVector[l1, skalartau - łP10P9start/(2 l1)];

tauP10st1st21plusstart =
łP10P9start/l1 + tauP93third[l1, rbp, xbp, ybp];
\[CurlyPhi]1drugi[t_] :=
tauP10st1st21plus[t] - tauP10st1st21plusstart + \[CurlyPhi]start2
\[CurlyPhi]w[t] = {0, 0, \[CurlyPhi]1drugi[t]};
skalar\[CurlyPhi]1drugi = versorz.\[CurlyPhi]w[t];
l2w = angleVector[l2, skalar\[CurlyPhi]1drugi];
l5 = l2 Cos[\[CurlyPhi]start2];
\[Beta]w[t_] := {0, 0, \[Beta][t]}
skalarbeta = versorz.\[Beta]w[t];
l4w = angleVector[l4, skalarbeta];
(*velocities of specific points:*)
lcbkw = (1/(m1 + m2)) (m1 l1w + m2 l2w)/
2 ;(* promien wektor Subscript[OC, 2]*)
vcbkw = Cross[D[tauw[t], t], lcbkw];
vP1w = D[P1[t], t];
lBC3w = l4w ((m4/2 + mD)/(m4 + mD));

vBC3w = Cross[D[\[Beta]w[t], t], lBC3w];

vBw = Cross[D[\[CurlyPhi]w[t], t], l2w];

vC3w = vBw + vBC3w;
(*heights of specific points:*)
versory = UnitVector[3, 2];
yC2 = Dot[versory, lcbkw];
ypw[t_] := Dot[versory, P1[t]]
lC3w = lBC3w + l2w;
yC3 = Dot[versory, lC3w];
(*elongation of spring*)
lambda[t_] := ypw[t] - ypw[t2]

(*moments of inertia*)
IIb = 1/3 m2 l2^2;
IIk = 1/3 m1 l1^2;

IIbk = IIb + IIk;
IIDwzgB = mD l4^2;

IIrzwzgB = (m4 l4^2)/3;
IIDrz = IIDwzgB + IIrzwzgB - (m4 + mD) (modul[lBC3w])^2;

(*potential energies with problematic Vspr component*)
VP1 = mpw g ypw[t];
Vspr = 1/2 kspr lambda[t]^2;
VC2 = (m1 + m2) g yC2;
VC3 = (m4 + mD) g yC3;
V = VP1 + Vspr + VC2 + VC3;
(*kinetic energies:*)
TP1r = 0;
TC2r = (IIbk modul[D[tauw[t], t] ]^2)/2;
TC3r = (IIDrz modul[D[\[Beta]w[t], t] ]^2)/2;
TP1p = (mpw modul[vP1w]^2)/2;
TC2p = ((m1 + m2) modul[vcbkw]^2)/2;
TC3p = ((m4 + mD) modul[vC3w]^2)/2;
T = TP1p + TC2p + TC3p + TP1r + TC2r + TC3r;
(*Lagrangian equations *)
Q = 0;
genCoords = {tauP10st1st21plus[t], \[Beta][t]};
LagrangianEquations[T_, V_, Q_: 0, genCoords_List] :=
Module[{L = T - V}, (D[D[L, D[#, t]], t] == Q + D[L, #]) & /@
genCoords]
eqLagrangeIIwSprgphaseIIst21 = LagrangianEquations[T, V, Q, genCoords];
(*data*)
l1 = 0.6096;
l2 = 2.4384;
rbp = 0.06;
xbp = 0.7096;
ybp = -0.1;
CC = 0;
ls0 = 0.5;
lsmax = 0.625;
ll2 = 0.4156500000000001;
rb2 = 0.12;
ybpP1tstart = 0.2;
t1 = 0;
\[CurlyPhi]start = 0.785398;
\[CurlyPhi]start2 = 5.497787;
rb5 = 0.3;
rb6 = 0.4;
rb7 = 0.1;
rb4 = 0.1;
yb4 = 0.5;
\[Beta]4 = (-3*Pi)/4;
\[Omega]\[Beta]3 = -0.5;
l4 = 2.1336;
l5 = 1.72421;
mb = 24.5455;
m1 = 4.9091000000000005;
m2 = 19.636400000000002;
mD = 3.63636;
m4 = 0;
g = 9.81;
mpw = 363.636;
time = 0.6;
kspr = 1000;
t2 = 0.5717686223842463;
(*boundary conditions*)
intercphaseIIst1st21tau = {tauP10st1st21plus[t2] ==
2.1275269499525282,
Derivative[1][tauP10st1st21plus][
t2] == -4.645459820712312, \[Beta][t2] == -0.5484939695351693,
Derivative[1][\[Beta]][t2] == -9.20429489396812};
(*solution*)
ndfazaIIst21Spr =
NDSolve[{eqLagrangeIIwSprgphaseIIst21,
intercphaseIIst1st21tau}, {tauP10st1st21plus, \[Beta]
}, {t, t2, 1}, Method -> {"EquationSimplification" -> "Residual"}]


Nice to have the ODEs and ICs separate, because it makes separate analyses easier:

Cases[eqLagrangeIIwSprgphaseIIst21,
Derivative[_][_][_] | tauP10st1st21plus[_] | β[_],
Infinity] // DeleteDuplicates
(*
{tauP10st1st21plus[t],
Derivative[1][tauP10st1st21plus][t],
Derivative[2][tauP10st1st21plus][t],
\[Beta][t],
Derivative[1][\[Beta]][t],
Derivative[2][\[Beta]][t],
tauP10st1st21plus[0.571769]}  <-- N.B.
*)

Cases[{intercphaseIIst1st21tau},
Derivative[_][_][_] | tauP10st1st21plus[_] | β[_],
Infinity] // DeleteDuplicates
(*
{tauP10st1st21plus[0.571769],
Derivative[1][tauP10st1st21plus][0.571769],
β[0.571769],
Derivative[1][β][0.571769]}
*)


We can see that the initial condition appears (partly) in the ODEs (first output above). InternalProcessEquationsSeparateEquations separates equations into differential, algebraic, and constraint (e.g. BCs) equations. So probably the appearance of tauP10st1st21plus[t2] in the ODEs causes the problem. It's easy to solve the ICs and substitute the solutions into the ODEs. Luckily, we can also solve the system for the leading derivatives so that we can avoid "EquationSimplification"->"Residual":

ndfazaIIst21Spr = NDSolve[{
Solve[eqLagrangeIIwSprgphaseIIst21,
{tauP10st1st21plus''[t], β''[t]}] /.
Rule -> Equal /.
First@Solve[intercphaseIIst1st21tau,
{tauP10st1st21plus[t2],
tauP10st1st21plus'[t2], β[t2], β'[t2]}]
, intercphaseIIst1st21tau}, {tauP10st1st21plus, β}
, {t, t2, 1}
(*,Method->{"EquationSimplification"->"Residual"}*)];

ndfazaIIst21Spr // Values // Flatten // ListLinePlot


• Hello, thank you for your explanation. I havent understood all of it in deep so far however I can see that your solution is wrong - I can obtain the same result when I omit that troublesome Vspr component. Jun 8, 2022 at 16:24
• @lodzki I can't comment on the in/correctness of your model. The following variant gives the same solution, which suggests the solution to the model as given in your code is accurately computed both above and as follows: NDSolve[{eqLagrangeIIwSprgphaseIIst21 /. First@Solve[intercphaseIIst1st21tau, {tauP10st1st21plus[t2], tauP10st1st21plus'[t2], \[Beta][t2], \[Beta]'[t2]}], intercphaseIIst1st21tau}, {tauP10st1st21plus, \[Beta]}, {t, t2, 1}, Method -> {"EquationSimplification" -> "Residual"}]. Jun 8, 2022 at 16:48
• The solution is also stable under increasing WorkingPrecision. Maybe there is a mistake in the setup of the ODE system? Jun 8, 2022 at 16:50
• I hope you're right - Probably I should set another generalized coordinate for spring - Vspr is for spring energy. Jun 8, 2022 at 18:35
• I will set your post as an answer because it is a solution to the problem I stated in my question (the DDE case) nevertheless NDSolve still computes Vspr wrong- it looks like inner flaw of the solver. Now I know how to omit that obstacle but still dont why such a thing occurs - I will start a new thread about it. Jun 15, 2022 at 11:33