Improvement of code precision using NDSolve for Differential-Algebraic equation

Question

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.

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

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

The remaining dependent variables are identically zero.