# Evaluation Control for Equation Input to NDSolve

I am trying to numerically integrate a system of equations using NDSolve and am having issues with symbolic matrix inversion taking a long time. The actual system is much more complicated, but I have made a simple example below. Let's assume the system can be written in the form:

q''[t] = Inverse[a[q[t],q'[t]].q'[t]


where q[t] is a 4x1 vector of variables, and a is a function of the state variables. If we specify q' and q we should be able to solve for q[t]. The problem is a[q[t],q'[t]] is a very large symbolic 4x4 matrix so trying to solve the inversion symbolically is very time consuming. Is there a way to force the matrix inversion to hold off (using something like Evaluation Control) until the symbolic state variables are substituted with the numerical values by NDSolve? It looks like a related question was asked before here but without a good response.

Example source code below:

n = 3 (* Sets exponent for the "a" matrix, fast if n=1, slow for n >=3*)

q = {{q1[t]}, {q2[t]}, {q3[t]}, {q4[t]}} (* configuration variables*)
dq = D[q, t] (* velocity variables*)
ddq = D[dq, t] (* acceleration variables*)

(* Setting up initial conditions*)
q0 = Thread[Flatten@q == 0.1] /. t -> 0
dq0 = Thread[Flatten@dq == -0.1] /. t -> 0

(* Large matrix that needs to be inverted A = Inverse[a[q,dq]]*)
A[q_, dq_, n_] := Inverse[
{{2*Cos[dq[]], 2*Sin[q[]], q[] + Sin[q[]], Cos[dq[]] + Sin[dq[]]},
{2*Sin[q[]], 2*q[], Cos[dq[]] + Sin[q[]], Cos[dq[]] + dq[]},
{q[] + Sin[q[]], Cos[dq[]] + Sin[q[]], q[] + Sin[q[]], Cos[dq[]] + q[]},
{Cos[dq[]] + dq[], Cos[dq[]] + dq[], Cos[dq[]] + q[], dq[] + Tan[dq[]]}}^n]

{tvar, RHS} = Timing[Flatten[A[Flatten@q, Flatten@dq, n].dq]];
Print[tvar] (* print execution time for evaluating A matrix*)

ddqEqs = Thread[Flatten@ddq == RHS]; (* Setting up equations*)
NDSolve[Join[ddqEqs, q0, dq0], Flatten@q, {t, 0, 1}] (* Run numerical integration*)


I recommend that you vectorize your ODE, and use LinearSolve instead of Inverse. To do so, let's first rewrite your matrix using Indexed:

m = Quiet @ {
{2*Cos[dq[]],2*Sin[q[]],q[]+Sin[q[]],Cos[dq[]]+Sin[dq[]]},
{2*Sin[q[]],2*q[],Cos[dq[]]+Sin[q[]],Cos[dq[]]+dq[]},
{q[]+Sin[q[]],Cos[dq[]]+Sin[q[]],q[]+Sin[q[]],Cos[dq[]]+q[]},
{Cos[dq[]]+dq[],Cos[dq[]]+dq[],Cos[dq[]]+q[],dq[]+Tan[dq[]]}
} /. Part->Indexed;


Since Indexed is not supported by TeXForm, I will replace Indexed with Subscript to show it in MathJax:

m /. Indexed -> Subscript //TeXForm


$$\left( \begin{array}{cccc} 2 \cos \left(\text{dq}_{\{1\}}\right) & 2 \sin \left(q_{\{4\}}\right) & q_{\{3\}}+\sin \left(q_{\{4\}}\right) & \sin \left(\text{dq}_{\{4\}}\right)+\cos \left(\text{dq}_{\{1\}}\right) \\ 2 \sin \left(q_{\{4\}}\right) & 2 q_{\{1\}} & \cos \left(\text{dq}_{\{1\}}\right)+\sin \left(q_{\{3\}}\right) & \text{dq}_{\{1\}}+\cos \left(\text{dq}_{\{1\}}\right) \\ q_{\{3\}}+\sin \left(q_{\{4\}}\right) & \cos \left(\text{dq}_{\{1\}}\right)+\sin \left(q_{\{3\}}\right) & q_{\{4\}}+\sin \left(q_{\{4\}}\right) & \cos \left(\text{dq}_{\{1\}}\right)+q_{\{2\}} \\ \text{dq}_{\{4\}}+\cos \left(\text{dq}_{\{1\}}\right) & \text{dq}_{\{1\}}+\cos \left(\text{dq}_{\{1\}}\right) & \cos \left(\text{dq}_{\{1\}}\right)+q_{\{2\}} & \text{dq}_{\{4\}}+\tan \left(\text{dq}_{\{4\}}\right) \\ \end{array} \right)$$

Now, we can define the RHS of the equation (the somewhat convoluted syntax below is to fix scoping issues):

rhs[q_List, dq_List] = With[{m=m}, Unevaluated @ LinearSolve[m, dq]];


Because of the List pattern, rhs will only evaluate when fed list arguments:

rhs[q, dq]
rhs[N@{1, 2, 3, 4}, N@{4, 3, 2, 1}]


rhs[q, dq]

{-0.906334, -0.657835, 0.728103, 0.99123}

Now, we are ready to use NDSolve (or rather NDSolveValue since it is more convenient):

sol = NDSolveValue[
{
q''[t] == rhs[q[t], q'[t]],
q == ConstantArray[.1, 4],
q' == ConstantArray[-.1, 4]
},
q,
{t, 0, 1}
];


Visualization:

Plot[sol[t], {t, 0, 1}] • Thanks for the helpful response. I think this partially solves what I'm trying to do, I'll work through it and report back with what I figure out. – Zack May 24 '19 at 19:44