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.
My code is below. Please help me.
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)
(*radius-vectors of specific points:*)
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"}]
