6
$\begingroup$

I am trying to find a numerical solution to a set of coupled ODEs with NDSolve. Let us say $\boldsymbol{X}$ is a vector, and $\boldsymbol{F}$ is a non-linear vector function of $\boldsymbol{X}$: $$\dot{\boldsymbol{X}}=\boldsymbol{F}(\boldsymbol{X}),\qquad\dot{\boldsymbol{X}}(0)=\boldsymbol{X}_{0}$$.

In a nutshell, my problem is this: $\boldsymbol{F}$ is basically infeasible to evaluate analytically (if you pass it a symbolic list $\left\{ X_{1},X_{2},\ldots,X_{n}\right\} $ as an argument). However, if you pass it a list of numbers, $\boldsymbol{F}$ is no problem to evaluate. NDSolve seems to try to evaluate things analytically (I think), preventing it from calculating a result in feasible time. I'm looking for a way to tell NDSolve to stop trying any analytical evaluation of $\boldsymbol{F}$.

Here is an example just to illustrate the point. Let's say $\boldsymbol{F}$ involves some nasty matrix inversions: $$\boldsymbol{F}(\boldsymbol{X})=\left(\boldsymbol{1}+\left(\boldsymbol{1}+\boldsymbol{X}\otimes\boldsymbol{X}\right)^{-1}\right)^{-1}\boldsymbol{X}$$

ndim = 10;
Xvec[t_] = Table[x[k][t], {k, 1, ndim}];
F[xvec_] := 
 Nest[Inverse[IdentityMatrix[ndim] + #] &, TensorProduct[xvec, xvec], 
   2].xvec

In ndim = 10 dimensions, if you pass a symbolic argument just evaluating $\boldsymbol{F}$ takes longer than I have patience for:

    AbsoluteTiming[F[Xvec[t]];]
(* takes too long! *)

If instead you pass a numerical argument, it takes no time at all to evaluate:

AbsoluteTiming[F[RandomReal[{-5, 5}, ndim]]]

(* OUTPUT: {0.0000794, {3.72758, -2.68105, -1.71283, -3.62467, -3.39943, \
-1.28462, 2.47934, 0.216207, -0.944914, 1.36631}} *)

I think this problem carries over to NDSolve. The following calculation does not finish in feasible time for ndim = 10:

Xvec0 = RandomReal[{-5, 5}, ndim]; (* initial conditions *)
tend = 10;
sol = NDSolve[{Thread[Xvec'[t] == F[Xvec[t]]], 
   Thread[Xvec[0] == Xvec0]}, Xvec[t], {t, 0, tend}]

(* takes too long! *)

Is there a way I can instruct NDSolve to suppress analytical calculation and treat $\boldsymbol{F}$ with numerical arguments?

An unsuccessful solution attempt with NumericQ

I tried defining F with a pattern that checks for vector arguments (VectorQ) and that the argument is numerical (NumericQ):

ndim = 10;
Xvec[t_] = Table[x[k][t], {k, 1, ndim}];
F[xvec_?(VectorQ[#, NumericQ] &)] := 
 Module[{}, 
  Nest[Inverse[IdentityMatrix[ndim] + #] &, 
    TensorProduct[xvec, xvec], 2].xvec]

The behaviour is as expected: Now, symbolic arguments no longer evaluate

AbsoluteTiming[F[Xvec[t]]]

(* OUTPUT: {0.0000106, 
 F[{x[1][t], x[2][t], x[3][t], x[4][t], x[5][t], x[6][t], x[7][t], 
   x[8][t], x[9][t], x[10][t]}]} *)

... and numerical arguments evaluate quickly, just as previously.

However, NDSolve does not evaluate F, sadly:

Xvec0 = RandomReal[{-5, 5}, ndim]; (* initial conditions *)
tend = 10;
sol = NDSolve[{Thread[Xvec'[t] == F[Xvec[t]]], 
   Thread[Xvec[0] == Xvec0]}, Xvec[t], {t, 0, tend}]

(* OUTPUT: NDSolve[{{Derivative[1][x[1]][t] == F[{x[1][t], x[2][t], x[3][t], ... *)

I am hoping for some suggestions to make NDSolve calculate a solution in feasible time.

$\endgroup$
  • 2
    $\begingroup$ You're almost there, just remove Xvec[t_] = Table[x[k][t], {k, 1, ndim}]; and Threads. $\endgroup$ – xzczd Aug 5 at 12:37
8
$\begingroup$

It is strange that NDSolve does not produce warning messages but the reason behind the failure is that F[{x[1][t], ...}] does not evaluate to an explicit List when supplied a symbolic argument. So, when you evaluate

Thread[Xvec'[t] == F[Xvec[t]]]

you end up with all x[k]'[t] equated to the same F[{x[1][t], ...}], each of which will later evaluate to a vector when supplied a numeric argument. This makes NDSolve confused because you have scalars in the initial conditions for x[k][0].

You don't have to use explicit vectors at all. NDSolve can treat an array as a single point:

ClearAll[ndim, Xvec0, F, x];
ndim = 10;
Xvec0 = RandomReal[{-5, 5}, ndim];
F[xvec_?(VectorQ[#, NumericQ] &)] := 
Module[{}, 
  Nest[Inverse[IdentityMatrix[ndim] + #] &, TensorProduct[xvec, xvec], 2] . xvec]
NDSolve[{x'[t] == F[x[t]], x[0] == Xvec0}, x, {t, 0, 1}]
$\endgroup$
  • $\begingroup$ Thanks, this works, and your explanation is very helpful (especially that F does not evaluate to a list with symbolic arguments). $\endgroup$ – Alexander Erlich Aug 5 at 15:11

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.