2

I wish to solve an ordinary differential equation of matrices. This is a minimal working example I am working towards nonlinear equations. The accuracy of the solution is important for me. I have hit two problems.

  1. Why does Mathematica have to set up differential-algebraic equations?
  2. For accuracy how can I use something like the options; InterpolationOrder->All, Method->"ImplicitRungeKutta"?

I generate a set of three differential equations, their initial conditions and then go on to solve them using NDSolve

ClearAll[d, v, t];
ff0 = {10, 15, 22};
mm0 = {1, 2, 2.5};
dr = {0.001, 0.01, 0.015};
nn = Length[ff0];
SeedRandom[123];
R0 = Orthogonalize[RandomReal[{-1, 1}, {nn, nn}]];
R = Orthogonalize[RandomReal[{-1, 1}, {nn, nn}]];
cc = R0 . DiagonalMatrix[2 dr 2 \[Pi] ff0 ] . Transpose[R0];
mm = R . DiagonalMatrix[mm0] . Transpose[R];
kk = R . DiagonalMatrix[(2 \[Pi] ff0)^2 mm0] . Transpose[R];

(* Initial conditions *) dic = ConstantArray[0, nn]; dic[[1]] = 1; vic = ConstantArray[0, nn];

(* Solve differential equations *) te = 30; {disp, vel} = {d, v} /. First@NDSolve [{ d'[t] == v[t], mm . v'[t] + cc . v[t] + kk . d[t] == 0, d[0] == dic, v[0] == vic}, {d, v}, {t, 0, te}]

This does work and gives me answers but the following message appears

NDSolve::ntdvdae: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations.

Why is Mathematica doing this?

If I try for high accuracy using

(* Solve differential equations with high accuracy *)
te = 30;
{disp, vel} = {d, v} /. First@NDSolve [{
     d'[t] == v[t],
     mm . v'[t] + cc . v[t] + kk . d[t] == 0,
     d[0] == dic, v[0] == vic},
    {d, v},
    {t, 0, te},
    InterpolationOrder -> All,
    Method -> "ImplicitRungeKutta"]

Then I get a string of messages with these being key

NDSolve::ntdvdae: Cannot solve to find an explicit formula for the derivatives. NDSolve will try solving the system as differential-algebraic equations.

NDSolve::nodae: The method NDSolve`ImplicitRungeKutta is not currently implemented to solve differential-algebraic equations. Use Method -> Automatic instead.

I am using version 14.0.0 for windows.

Thanks for any help

Hugh
  • 16,387
  • 3
  • 31
  • 83
  • 1
    Use mm1 = Inverse[mm]; sol = NDSolve[{d'[t] == v[t], v'[t] + mm1 . (cc . v[t]) + mm1 . (kk . d[t]) == 0, d[0] == dic, v[0] == vic}, {d, v}, {t, 0, te}, Method -> "ImplicitRungeKutta"] – Alex Trounev Jan 28 '24 at 00:09
  • Seems that the pre-processor of NDSolve isn't strong enough. This system can be solved with MassMatrix method to avoid calculating Inverse[mm] in principle (MassMatrix sometimes works better than Residual, and it works with some of the ODE solvers: https://mathematica.stackexchange.com/a/277822/1871 ), but NDSolve just can't do it when v and d are implicit vectors. The only solution I can think out at the moment is to use explicit vector to represent v and d, the pre-processing will be slower then, but at least we can now adjust the ODE solver freely.
  • – xzczd Jan 28 '24 at 04:43
  • 2
  • Are you sure you need ImplicitRungeKutta? According to my limited experience, setting a proper MaxStepSize is usually quite enough to achieve good accuracy.
  • – xzczd Jan 28 '24 at 04:44
  • @AlexTrounev Thanks that does work. I will have to think if this is applicable when I go forward to nonlinear systems. – Hugh Jan 28 '24 at 21:26
  • @Hugh It could be better to show your nonlinear system. – Alex Trounev Jan 28 '24 at 21:37
  • @xzczd I will investigate controlling the step size. I was following the very successful output from this post – Hugh Jan 28 '24 at 21:41