I'm trying to solve a system of 24 non-linear Differential-Algebraic equations (DAE). I'm using the command NDSolve in Mathematica to solve this system, using this command, the error is too much large. I want to improve the precision of the code, for this I was trying different methods in NDSolve command. But, Mathematica is unable to solve. I'm getting the error:

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

I want to use the Implicit-Runge-Kutta method or projection method to improve my results.

If I used these methods in a system of ODE's in NDSolve command, mathematica is able to give output.

Just as an example to test the code, I'm posting here some short example:

NDSolve[{x'[t] == -y[t], y'[t] == x[t], x[0] == 0.1, y[0] == 0}, {x, 
  y}, {t, 0, 100}, 
 Method -> {"FixedStep", 
   Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 10, 
     "ImplicitSolver" -> {"Newton", AccuracyGoal -> MachinePrecision, 
       PrecisionGoal -> MachinePrecision, 
       "IterationSafetyFactor" -> 1}}}, StartingStepSize -> 1/10]

I'm able to obtain the output of the above system using Implicit-Runge-Kutta method, but If I use DAE system, I'm not able to get output, for example:

NDSolve[{x'[t] - y[t] == Sin[t], x[t] + y[t] == 1, x[0] == 0}, {x, 
  y}, {t, 0, 10}, 
 Method -> {"FixedStep", 
   Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 10, 
     "ImplicitSolver" -> {"Newton", AccuracyGoal -> 15, 
       PrecisionGoal -> 50, "IterationSafetyFactor" -> 1}}}, 
 StartingStepSize -> 1/10]

Can anyone help me please, how can I solve such a DAE system with NDSolve command using some implicit method, like Implicit-Runge-Kutta method?

Should I convert this DAE system into ODE's, if yes then how can we convert such a system into a system of ordinary differential equations?

Actually, I'm working in General Relativity, here to apply the method as for the above example is not simple. I'm still not able to solve the system. I'm posting here my system of DAE equations.

  • $\begingroup$ Comments are not for extended discussion; this conversation has been moved to chat. $\endgroup$
    – Kuba
    Oct 13, 2020 at 7:28

1 Answer 1


MichaelE2 already has answered the question in a comment: To use Method -> "ImplicitRungeKutta", differentiate the second equation and add a corresponding boundary condition for y. However, the OP expressed the concern that doing so might produce an inaccurate answer. Out of curiosity, I tried it. So, the following actually is an extended comment.

It is easy to determine the accuracy of any numerical solution to the system of equations, because a symbolic solution exists.

sa = DSolveValue[{x'[t] - y[t] == Sin[t], x[t] + y[t] == 1, x[0] == 0}, 
    {x[t], y[t]}, {t, 0, 10}];
(* {1/2 (2 - E^-t - Cos[t] + Sin[t]), 1/2 (E^-t + Cos[t] - Sin[t])} *)

Then, applying the approach recommended by MichaelE2,

NDSolveValue[{x'[t] - y[t] == Sin[t], x'[t] + y'[t] == 0, x[0] == 0, y[0] == 1}, 
   {x[t], y[t]}, {t, 0, 10}, Method -> "ImplicitRungeKutta", InterpolationOrder -> All];

yields excellent accuracy.

Plot[Evaluate[% - sa], {t, 0, 10}, PlotRange -> All, ImageSize -> Large, 
    AxesLabel -> {t, "x,y"}, LabelStyle -> {15, Bold, Black}]

enter image description here

Note that InterpolationOrder -> All is needed to eliminate spurious oscillations in the InterpolationFunction of order 10^-5. Whether this approach can be used in the 24-equation system mentioned by the OP depends on the details of those equations, which I have requested.

Incidentally, I find it surprising that NDSolve does not simplify the original DAE system to eliminate y[t] and numerically integrate the resulting ODE in x[t], instead of terminating when Method -> "ImplicitRungeKutta" is employed.

Addendum: Solution to set of 24 nonlinear equations

NDSolve misinterprets the system of enormous equations recently added to the question as a DAE system due to

Vup = {Vt[τ], Vr[τ], Vθ[τ], Vϕ[τ]};

These four quantities are, in fact, simply names for expressions and should be renamed as

Vup = {Vt, Vr, Vθ, Vϕ};

The code giving them values then becomes

{Vt, Vr, Vθ, Vϕ} = NN (Wvec + Pup) /. t -> t[τ] /. r -> r[τ] /. θ -> θ[τ] /. ϕ -> ϕ[τ];

instead of the expression for EQ4. Of course, EQ4 then must be deleted from the subsequent expression for EOM. The code leading to EOM also has an error somewhere, which I corrected rather inelegantly by inserting after the expression for EOM the further line of code,

EOM = EOM /. z_[τ][τ] -> z[τ];

With these changes NDSolve successfully runs until r[τ] decreases to 2, the event horizon. Specifically,

var = Through[{t, r, θ, ϕ, Pt, Pr, Pθ, Pϕ, Stt, Str, Stθ, Stϕ, 
    Srt, Srr, Srθ, Srϕ, Sθt, Sθr, Sθθ, Sθϕ, Sϕt, Sϕr, Sϕθ, Sϕϕ}[τ]];
NDSolveValue[Flatten[Join[{EOM, INT}]], var, {τ, 0, 1000}, 
    Method -> {"ImplicitRungeKutta"}];

terminates with NDSolveValue::ndsz at τ = 37.771696. A plot of the first eight variables then is,

Plot[Evaluate[%[[;; 8]]], {τ, 0, 37.77169}, PlotRange -> {Automatic, 8}, 
    ImageSize -> Large, PlotLegends -> Placed[ToString /@ var, {.35, .6}]]

enter image description here

The remaining dependent variables are identically zero.

  • 1
    $\begingroup$ Er… what do you mean by "I find it surprising that NDSolve does not simplify the system before starting the numerical integration"? $\endgroup$
    – xzczd
    Oct 12, 2020 at 5:48
  • 3
    $\begingroup$ @xzczd Good question. I have clarified my final sentence. $\endgroup$
    – bbgodfrey
    Oct 12, 2020 at 13:17
  • $\begingroup$ @bbgodfrey, thank for your great efforts, how can be the dependent variables zero? Even the posted code was giving the solution for all 24 variables. I used the approach you have mentioned but Mathematica is giving error:: DSolve::underdet: here are more dependent variables, than equations, so the system is underdetermined. I think we should adopt some different command together with the ImplicitRungeKutta command to solve such a system. I was following this, link, but unfortunately was not successful. $\endgroup$
    – MMS
    Oct 13, 2020 at 5:21
  • $\begingroup$ @bbgodfrey thanks a lot for your great efforts, now I'm able to run the code without any error. Actually, I changed the initial conditions, now I am able to obtain the solutions for all 24 variables. $\endgroup$
    – MMS
    Oct 13, 2020 at 19:46
  • 1
    $\begingroup$ @MMS You are welcome. I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. $\endgroup$
    – bbgodfrey
    Oct 13, 2020 at 20:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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