4
$\begingroup$

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"}]
Improve this question
$\endgroup$

1 Answer 1

Reset to default
4
$\begingroup$

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). Internal`ProcessEquations`SeparateEquations 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

enter image description here

Improve this answer
$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$ Commented Jun 8, 2022 at 16:24
  • $\begingroup$ @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"}]. $\endgroup$ Commented Jun 8, 2022 at 16:48
  • $\begingroup$ The solution is also stable under increasing WorkingPrecision. Maybe there is a mistake in the setup of the ODE system? $\endgroup$ Commented Jun 8, 2022 at 16:50
  • $\begingroup$ I hope you're right - Probably I should set another generalized coordinate for spring - Vspr is for spring energy. $\endgroup$ Commented Jun 8, 2022 at 18:35
  • $\begingroup$ 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. $\endgroup$ Commented Jun 15, 2022 at 11:33

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.