# How to update dependent variable at each time steps of NDSolve

I want to numerically solve the following ODEs in Mathematica.

$\ddot{x}(t)+y(t)=\cos t$ and $\ddot{y}(t)+x(t)=\mathrm{e}^t$, subjected to a constraint equation $C(x,y)=x^2+y=0$.

Initial conditions are $x(0)=1$ and $\dot{x}(0)=0$.

The objective is to correct $x(t)$ and $y(t)$ at each time step according to following equation

$\mathbf{q}^\mathrm{c}=\mathbf{q}^\mathrm{u}-C_\mathbf{q}^\mathrm{T}(C_\mathbf{q}C_\mathbf{q}^\mathrm{T})^{-1}C(\mathbf{q}^\mathrm{u})$ which we have obtained from the constraint equation. Here, $\mathbf{q}=\{x(t),y(t)\}^\mathrm{T}$. $\mathbf{q}^\mathrm{c}$ denotes the corrected coordinates while $\mathbf{q}^\mathrm{u}$ represents the uncorrected ones obtained by solving the ODEs at a particular time step. $C_\mathbf{q}=\{\partial C/\partial x\quad\partial C/\partial y\}$.

After correcting $\mathbf{q}$, we use the same (i.e., $\mathbf{q}^\mathrm{c}$) to correct the first derivatives. Therefore, we use the following equation $\mathbf{v}^\mathrm{c}=\mathbf{v}^\mathrm{u}-C_\mathbf{q}^\mathrm{T}(C_\mathbf{q}C_\mathbf{q}^\mathrm{T})^{-1}\dot{C}(\mathbf{q}^\mathrm{c},\mathbf{v}^\mathrm{u})$. Here, $\mathbf{v}=\{\dot{x}(t),\dot{y}(t)\}^\mathrm{T}$ and $\dot{C}=2x\dot{x}+\dot{y}$.

I am not familiar with updating variables while using NDSolve at every time steps.

The following is the code for dynamic simulation of a slider crank mechanism, on which I want to apply the correction method. Earlier I had posted a simple mathematics which was, unfortunately, not very clear and also created a bit of confusion.

The dynamics of a multibody system can be solved by various methods. One of them is to the write the equations of motion (i.e., the ODEs) by increasing the degrees of freedom (or the coordinates) of the multibody system. The constraints are incorporated by introducing Lagrange multipliers. This method reduces the computational effort by a great amount but at the cost of constraint violation. With the passage of time, the errors in determining the coordinates (i.e. $\mathbf{q}$) which come from the numerical integration of ODEs accumulate and after a certain period of time the coordinates fail to satisfy the constraint equations. To handle the issue various correction methods are applied, one of which has been provided in my question. This corrections are done at each time steps of the numerical integration. Below I have provided results of $\theta_3$ and $\lvert \mathbf{R}_2 \rvert$ only to explain the phenomenon of constraint violation.

It can be seen that after a certain period of time the response of $\theta_3$ starts deviating from the equilibrium line (orange one). And also $\lvert \mathbf{R}_2 \rvert$ should remain constant throughout the simulation, which did not happen in practical.

 ClearAll["Global*"]

BlockDiagonalMatrix[b:{__?MatrixQ}] := Module[{r, c, n = Length[b], i, j}, {r, c} = Transpose[Dimensions /@ b]; ArrayFlatten[Table[If[i == j, b[[i]], ConstantArray[0, {r[[i]], c[[j]]}]], {i, n}, {j, n}]]]
(*-------------------------------Generalized Coordinates------------------------------------------------------------------------------------------------------------------------------------*)
R2[t] = Transpose[{{Rx2[t], Ry2[t]}}];
R3[t] = Transpose[{{Rx3[t], Ry3[t]}}];
R4[t] = Transpose[{{Rx4[t], Ry4[t]}}];
q2 = ArrayFlatten[Transpose[{{R2[t], \[Theta]2[t]}}]];
q3 = ArrayFlatten[Transpose[{{R3[t], \[Theta]3[t]}}]];
q4 = ArrayFlatten[Transpose[{{R4[t], \[Theta]4[t]}}]];
q = ArrayFlatten[Transpose[{{q2, q3, q4}}]];
(*-----------------------Time Derivative of Generalized Coordinates--------------------------------------------------------------------------------------------------------------------------------------*)
dqdt = D[q, t];
(*-------------------------------------Data-------------------------------------------------------------------------------------------------------------------------------*)
m2 = 15; m3 = 30; m4 = 15;
l0 = 0.305; l2 = 1.67; l3 = 3.07; b = 0.15; a = 0.15/2; e = 1.58; g = 9.81; k = 8000; kO2 = 8380; kA3 = 6200; c = 100; cO = 0; cA = 0;
cB = 0;
J2 = (1/12)*m2*l2^2; J3 = (1/12)*m3*l3^2; J4 = (1/12)*m4*(a^2 + b^2);
kO = kO2*(l2^2/4)*Cos[\[Theta]2[t]]^2;
kA = kA3*(l3^2/4)*Cos[\[Theta]3[t] - \[Theta]2[t]]^2;
kB = 0;
(*--------------------------------Formation of Mass Matrix------------------------------------------------------------------------------------------------------------------------------------*)
M2 = DiagonalMatrix[{m2, m2, J2}];
M3 = DiagonalMatrix[{m3, m3, J3}];
M4 = DiagonalMatrix[{m4, m4, J4}];
M = BlockDiagonalMatrix[{M2, M3, M4}];
(*------------------------------Formation of Constraint Equations--------------------------------------------------------------------------------------------------------------------------------------*)
xBarO2 = -(l2/2); yBarO2 = 0; xBarA2 = l2/2; yBarA2 = 0; xBarA3 = -(l3/2); yBarA3 = 0; xBarB3 = l3/2; yBarB3 = 0;
xBarB4 = 0;
yBarB4 = 0;
uBarO2 = {{xBarO2}, {yBarO2}};
A2 = {{Cos[\[Theta]2[t]], -Sin[\[Theta]2[t]]}, {Sin[\[Theta]2[t]], Cos[\[Theta]2[t]]}};
A\[Theta]2 = D[A2, \[Theta]2[t]];
uBarA2 = {{xBarA2}, {yBarA2}};
uBarA3 = {{xBarA3}, {yBarA3}};
A3 = {{Cos[\[Theta]3[t]], -Sin[\[Theta]3[t]]}, {Sin[\[Theta]3[t]], Cos[\[Theta]3[t]]}};
A\[Theta]3 = D[A3, \[Theta]3[t]];
uBarB3 = {{xBarB3}, {yBarB3}};
uBarB4 = {{xBarB4}, {yBarB4}};
A4 = {{Cos[\[Theta]4[t]], -Sin[\[Theta]4[t]]}, {Sin[\[Theta]4[t]], Cos[\[Theta]4[t]]}};
A\[Theta]4 = D[A4, \[Theta]4[t]];
C1 = R2[t] + A2 . uBarO2;
C2 = R2[t] + A2 . uBarA2 - R3[t] - A3 . uBarA3;
C3 = R3[t] + A3 . uBarB3 - R4[t] - A4 . uBarB4;
C4 = {{Rx4[t] - e}, {\[Theta]4[t]}};
CE = ArrayFlatten[{{C1}, {C2}, {C3}, {C4}}];
(*---------------------------Formation of Cq Matrix-----------------------------------------------------------------------------------------------------------------------------------------*)
CqO = ArrayFlatten[{{IdentityMatrix[2], A\[Theta]2 . uBarO2, ConstantArray[0, {2, 2}], ConstantArray[0, {2, 1}], ConstantArray[0, {2, 2}], ConstantArray[0, {2, 1}]}}];
CqA = ArrayFlatten[{{IdentityMatrix[2], A\[Theta]2 . uBarA2, -IdentityMatrix[2], -A\[Theta]3 . uBarA3, ConstantArray[0, {2, 2}], ConstantArray[0, {2, 1}]}}];
CqB = ArrayFlatten[{{ConstantArray[0, {2, 2}], ConstantArray[0, {2, 1}], IdentityMatrix[2], A\[Theta]3 . uBarB3, -IdentityMatrix[2], -A\[Theta]4 . uBarB4}}];
CqBSld = ArrayFlatten[{{ConstantArray[0, {2, 2}], ConstantArray[0, {2, 1}], ConstantArray[0, {2, 2}], ConstantArray[0, {2, 1}], {{1, 0, 0}, {0, 0, 1}}}}];
Cq = ArrayFlatten[{{CqO}, {CqA}, {CqB}, {CqBSld}}];
CqSwap = Cq;
CqSwap[[All,{8, 9}]] = CqSwap[[All,{9, 8}]];
Cqi = CqSwap[[All,9]];
Cqd = CqSwap[[1 ;; 8,1 ;; 8]];
Cdi = -Inverse[Cqd] . Cqi;
Qd = Simplify[-ArrayFlatten[{Table[D[Cq . dqdt, q[[i]]], {i, 1, 9}]}] . dqdt];
(*-----------------------------External Force Components--------------------------------------------------------------------------------------------------------------------------------------*)
xBarF2 = NA; yBarF2 = NA; Fx2 = 0; Fy2 = 0; xBarF3 = NA; yBarF3 = NA; Fx3 = 0; Fy3 = 0; xBarF4 = 0; yBarF4 = b/2; Fx4 = 0; Fy4 = 0; \[Tau]2AE = 0; HzVal = 0;
\[CapitalOmega] = 2*Pi*HzVal;
\[Tau]2Ext = \[Tau]2AE*Cos[\[CapitalOmega]*t];
uBarF2 = {{xBarF2}, {yBarF2}};
uF2 = ArrayFlatten[{{A2 . uBarF2}, {0}}];
F2 = {{Fx2}, {Fy2}, {0}};
\[Tau]2 = Cross[Transpose[uF2][[1]], Transpose[F2][[1]]] . {{0}, {0}, {1}};
Qe2 = Transpose[ArrayFlatten[{{{F2[[1 ;; 2,1]]}, \[Tau]2Ext + \[Tau]2[[1]] + ((kO + k*l0^2)*((\[Theta]1[t] - \[Theta]2[t]) - (\[Theta]1Start - \[Theta]2Start)) + (cO + c*l0^2)*(Derivative[1][\[Theta]1][t] - Derivative[1][\[Theta]2][t])) -
(kA*((\[Theta]2[t] - \[Theta]3[t]) - (\[Theta]2Start - \[Theta]3Start)) + cA*(Derivative[1][\[Theta]2][t] - Derivative[1][\[Theta]3][t]))}}] + {{0, (-m2)*g, 0}}] /. {\[Theta]1[t] -> 0, Derivative[1][\[Theta]1][t] -> 0, \[Theta]1Start -> 0};
uBarF3 = {{xBarF3}, {yBarF3}};
uF3 = ArrayFlatten[{{A3 . uBarF3}, {0}}];
F3 = {{Fx3}, {Fy3}, {0}};
\[Tau]3 = Cross[Transpose[uF3][[1]], Transpose[F3][[1]]] . {{0}, {0}, {1}};
Qe3 = Transpose[ArrayFlatten[{{{F3[[1 ;; 2,1]]}, \[Tau]3[[1]] + (kA*((\[Theta]2[t] - \[Theta]3[t]) - (\[Theta]2Start - \[Theta]3Start)) + cA*(Derivative[1][\[Theta]2][t] - Derivative[1][\[Theta]3][t]))}}] + {{0, (-m3)*g, 0}}];
uBarF4 = {{xBarF4}, {yBarF4}};
uF4 = ArrayFlatten[{{A4 . uBarF4}, {0}}];
F4 = {{Fx4}, {Fy4}, {0}};
\[Tau]4 = Cross[Transpose[uF4][[1]], Transpose[F4][[1]]] . {{0}, {0}, {1}};
Qe4 = Transpose[ArrayFlatten[{{{F4[[1 ;; 2,1]]}, \[Tau]4[[1]] + (kB*((\[Theta]3[t] - \[Theta]4[t]) - (\[Theta]3Start - \[Theta]4Start)) + cB*(Derivative[1][\[Theta]3][t] - Derivative[1][\[Theta]4][t]))}}] + {{0, (-m4)*g, 0}}];
Qe = ArrayFlatten[Transpose[{{Qe2, Qe3, Qe4}}]];
(*---------------------------Formation of Final Equations----------------------------------------------------------------------------------------------------------------------------------------*)
H\[Lambda]\[Lambda] = -Inverse[Cq . Inverse[M] . Transpose[Cq]];
Hq\[Lambda] = -Inverse[M] . Transpose[Cq] . H\[Lambda]\[Lambda];
Hqq = Inverse[M] + Inverse[M] . Transpose[Cq] . H\[Lambda]\[Lambda] . Cq . Inverse[M];
H\[Lambda]q = Transpose[Hq\[Lambda]];
RHS1 = Hqq . Qe + Hq\[Lambda] . Qd;
RHS2 = H\[Lambda]q . Qe + H\[Lambda]\[Lambda] . Qd;
(*-----------------------Setting Values of Coordinates of the Mechanism----------------------------------------------------------------------------------------------------------------------------------------*)
Ry4Start = 2.3;
DepVarStart = Solve[(Flatten[CE] /. {Ry4[t] -> Ry4Start}) == 0 && 0 < \[Theta]2[t] < Pi, {Rx2[t], Ry2[t], \[Theta]2[t], Rx3[t], Ry3[t], \[Theta]3[t], Rx4[t], \[Theta]4[t]}][[1]] /. {C[1] -> 0, C[2] -> 0};
R2Start = R2[t] /. DepVarStart;
\[Theta]2Start = \[Theta]2[t] /. DepVarStart;
R3Start = R3[t] /. DepVarStart;
\[Theta]3Start = \[Theta]3[t] /. DepVarStart;
R4Start = {{Rx4[t] /. DepVarStart}, {Ry4Start}};
\[Theta]4Start = \[Theta]4[t] /. DepVarStart;
Solve::ratnz
(*---------------------------------Initial Conditions-----------------------------------------------------------------------------------------------------------------------------------*)
Ry40 = 1.81478;
DepVarIni = Solve[(Flatten[CE] /. {Ry4[t] -> Ry40}) == 0 && 0 < \[Theta]2[t] < Pi, {Rx2[t], Ry2[t], \[Theta]2[t], Rx3[t], Ry3[t], \[Theta]3[t], Rx4[t], \[Theta]4[t]}][[1]] /. {C[1] -> 0, C[2] -> 0};
{{Rx20}, {Ry20}} = R2[t] /. DepVarIni;
\[Theta]20 = \[Theta]2[t] /. DepVarIni;
{{Rx30}, {Ry30}} = R3[t] /. DepVarIni;
\[Theta]30 = \[Theta]3[t] /. DepVarIni;
Rx40 = Rx4[t] /. DepVarIni;
\[Theta]40 = \[Theta]4[t] /. DepVarIni;
dRy4dt0 = 0;
dqddt0 = (Cdi /. Flatten[{Ry4[t] -> Ry40, DepVarIni}])*dRy4dt0;
{dRx2dt0, dRy2dt0, d\[Theta]2dt0, dRx3dt0, dRy3dt0, d\[Theta]3dt0, dRx4dt0, d\[Theta]4dt0} = dqddt0;
Solve::ratnz
(*-----------------------Development of ODEs--------------------------------------------------------------------------------------------------------------------------------------------*)
ode1 = Derivative[2][Rx2][t] == RHS1[[1,1]];
ode2 = Derivative[2][Ry2][t] == RHS1[[2,1]];
ode3 = Derivative[2][\[Theta]2][t] == RHS1[[3,1]];
ode4 = Derivative[2][Rx3][t] == RHS1[[4,1]];
ode5 = Derivative[2][Ry3][t] == RHS1[[5,1]];
ode6 = Derivative[2][\[Theta]3][t] == RHS1[[6,1]];
ode7 = Derivative[2][Rx4][t] == RHS1[[7,1]];
ode8 = Derivative[2][Ry4][t] == RHS1[[8,1]];
ode9 = Derivative[2][\[Theta]4][t] == RHS1[[9,1]];
intD1 = Derivative[1][Rx2][0] == dRx2dt0;
intD2 = Derivative[1][Ry2][0] == dRy2dt0;
intD3 = Derivative[1][\[Theta]2][0] == d\[Theta]2dt0;
intD4 = Derivative[1][Rx3][0] == dRx3dt0;
intD5 = Derivative[1][Ry3][0] == dRy3dt0;
intD6 = Derivative[1][\[Theta]3][0] == d\[Theta]3dt0;
intD7 = Derivative[1][Rx4][0] == dRx4dt0;
intD8 = Derivative[1][Ry4][0] == dRy4dt0;
intD9 = Derivative[1][\[Theta]4][0] == d\[Theta]4dt0;
int1 = Rx2[0] == Rx20;
int2 = Ry2[0] == Ry20;
int3 = \[Theta]2[0] == \[Theta]20;
int4 = Rx3[0] == Rx30;
int5 = Ry3[0] == Ry30;
int6 = \[Theta]3[0] == \[Theta]30;
int7 = Rx4[0] == Rx40;
int8 = Ry4[0] == Ry40;
int9 = \[Theta]4[0] == \[Theta]40;
system = {ode1, ode2, ode3, ode4, ode5, ode6, ode7, ode8, ode9, intD1, intD2, intD3, intD4, intD5, intD6, intD7, intD8, intD9, int1, int2, int3, int4, int5, int6, int7, int8, int9};
ti = 0; tf = 2000;
(*-------------------------Solutions-----------------------------------------------------------------------------------------------------------------------------------------*)
sol = NDSolve[system, {Rx2[t], Ry2[t], \[Theta]2[t], Rx3[t], Ry3[t], \[Theta]3[t], Rx4[t], Ry4[t], \[Theta]4[t]}, {t, ti, tf}][[1]];
(*--------------------Plotting the Results------------------------------------------------------------------------------------------------------------------------------------------------*)
Plot[Sqrt[Rx2[t]^2 + Ry2[t]^2] /. sol, {t, ti, tf}, PlotRange -> All, Frame -> True, FrameLabel -> {Automatic, "\!$$\*SuperscriptBox[\(R$$, $$2$$]\)"}, FrameStyle -> Directive[Black, Thick]]
Plot[Sqrt[Rx3[t]^2 + Ry3[t]^2] /. sol, {t, ti, tf}, PlotRange -> All, Frame -> True, FrameLabel -> {Automatic, "\!$$\*SuperscriptBox[\(R$$, $$3$$]\)"}, FrameStyle -> Directive[Black, Thick]]
Plot[{Rx4[t] /. sol, Ry4[t] /. sol}, {t, ti, tf}, Frame -> True, FrameLabel -> {Automatic, "\!$$\*\nStyleBox[SuperscriptBox[SubscriptBox[\"R\", \"x\"], \"4\"],\nFontWeight->\"Plain\"]$$\!$$\*\nStyleBox[\",\",\nFontWeight->\"Plain\"]$$\!$$\*\nStyleBox[\" \ \",\nFontWeight->\"Plain\"]$$\!$$\*\nStyleBox[SuperscriptBox[SubscriptBox[\"R\", \"y\"], \"4\"],\nFontWeight->\"Plain\"]$$"}, FrameStyle -> Directive[Black, Thick]]
Plot[{(180/Pi)*\[Theta]2[t] /. sol, 142.422}, {t, ti, tf}, Frame -> True, FrameLabel -> {Automatic, "\!$$\*SuperscriptBox[\(\[Theta]$$, $$2$$]\)"}, FrameStyle -> Directive[Black, Thick], PlotRange -> All]
Plot[{(180/Pi)*\[Theta]3[t] /. sol, 18.935}, {t, ti, tf}, Frame -> True, FrameLabel -> {Automatic, "\!$$\*SuperscriptBox[\(\[Theta]$$, $$3$$]\)"}, FrameStyle -> Directive[Black, Thick], PlotRange -> All]
Plot[{(180/Pi)*\[Theta]4[t] /. sol}, {t, ti, tf}, Frame -> True, FrameLabel -> {Automatic, "\!$$\*SuperscriptBox[\(\[Theta]$$, $$4$$]\)"}, FrameStyle -> Directive[Black, Thick], PlotRange -> All]


In the above code I have only used simple NDSolve to perform the integration (see Solutions) section. But I need to correct $\mathbf{q}$ at each time step according to the correction scheme mentioned earlier in my questions. I do not know how to incorporate this correction within NDSolve. $C$ is written as CE in 'Formation of Constraint Equation Section'. However, $\dot{C}$ is not given in the code. A simple example would also help. There is no need to spend musch time on this clumsy code.

• You say that you want to apply a correction based on a constraint equation. Could you not include the constraint as an additional equation and let NDSolve deal with it, rather than trying to do it by hand? What is your constraint equation like? Can you include some code in addition to the LaTeX expressions? – MarcoB Jun 26 '18 at 19:08
• I suspect your system is overdetermined: Substituting y[t]==-x[t]^2 into the two differential equations you will get two ode's in x[t] {x''[t] + y[t] == Cos[t], y''[t] + x[t] == Exp[t] } /. y -> (-x[#]^2 &) (*{-x[t]^2 + (x^\[Prime]\[Prime])[t] == Cos[t], x[t] - 2 Derivative[1][x][t]^2 - 2 x[t] (x^\[Prime]\[Prime])[t] == E^ t}*)` – Ulrich Neumann Jun 26 '18 at 20:00
• People here generally like users to post code as Mathematica code instead of just images or TeX, so they can copy-paste it. It makes it convenient for them and more likely you will get someone to help you. You may find this this meta Q&A helpful – Michael E2 Jun 26 '18 at 22:24