Evaluation Control for Equation Input to NDSolve

Question

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'[0] and q[0] 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[[1]]], 2*Sin[q[[4]]], q[[3]] + Sin[q[[4]]], Cos[dq[[1]]] + Sin[dq[[4]]]},
 {2*Sin[q[[4]]], 2*q[[1]], Cos[dq[[1]]] + Sin[q[[3]]], Cos[dq[[1]]] + dq[[1]]},
 {q[[3]] + Sin[q[[4]]], Cos[dq[[1]]] + Sin[q[[3]]], q[[4]] + Sin[q[[4]]], Cos[dq[[1]]] + q[[2]]}, 
 {Cos[dq[[1]]] + dq[[4]], Cos[dq[[1]]] + dq[[1]], Cos[dq[[1]]] + q[[2]], dq[[4]] + Tan[dq[[4]]]}}^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*)

Answer 1

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[[1]]],2*Sin[q[[4]]],q[[3]]+Sin[q[[4]]],Cos[dq[[1]]]+Sin[dq[[4]]]},
    {2*Sin[q[[4]]],2*q[[1]],Cos[dq[[1]]]+Sin[q[[3]]],Cos[dq[[1]]]+dq[[1]]},
    {q[[3]]+Sin[q[[4]]],Cos[dq[[1]]]+Sin[q[[3]]],q[[4]]+Sin[q[[4]]],Cos[dq[[1]]]+q[[2]]},
    {Cos[dq[[1]]]+dq[[4]],Cos[dq[[1]]]+dq[[1]],Cos[dq[[1]]]+q[[2]],dq[[4]]+Tan[dq[[4]]]}
} /. 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[0] == ConstantArray[.1, 4],
    q'[0] == ConstantArray[-.1, 4]
    },
    q,
    {t, 0, 1}
];

Visualization:

Plot[sol[t], {t, 0, 1}]