0
$\begingroup$

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*)
$\endgroup$
1
$\begingroup$

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}]

enter image description here

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Zack May 24 at 19:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.