Mathematica internal source code for a numerical method

Question

I read this (Six Reasons Why the Wolfram Language Is (Like) Open Source), which inspired me to re-investigate Mathematica's implementation of some numerical methods. I would like to see the code that evaluates, say, a differential equation via NDSolve[___,Method->((Some Method))].

I've tried the approaches given in the blog post (i.e., ResourceFunction["PrintDefinitions"][NDSolve]), but this only gives me a list of available NDSolve options. Some searching around revealed this beautiful package, Spelunking, but that seems to do more or less the same thing.

Spelunk[NDSolve]

{
 {{
   {Options[NDSolve] = {AccuracyGoal -> Automatic, 
      Compiled -> Automatic, DependentVariables -> Automatic, 
      DiscreteVariables -> {}, EvaluationMonitor -> None, 
      InitialSeeding -> {}, InterpolationOrder -> Automatic, 
      MaxStepFraction -> 1/10, MaxSteps -> Automatic, 
      MaxStepSize -> Automatic, Method -> Automatic, 
      NormFunction -> Automatic, PrecisionGoal -> Automatic, 
      StartingStepSize -> Automatic, StepMonitor -> None, 
      WorkingPrecision -> MachinePrecision}}
  }}
}

I've also attempted defining functions that call NDSolve with a specific method, but that did not reveal anything useful. I also found other Q/A pairs that pointed to Trace[] and TracePrint[], but those also didn't reveal anything substantive.

If Mathematica is (like) open source, how can I find the source for this example (or ones like it)?

Update: I found the option Trace[NDSolve[{x'[t] == t, x[0] == 0}, x, {t, 0, 6}], TraceInternal -> True] that opens up a can of back-end worms. This might be the best we can do (h/t this post).

Answer 1

I feel I've gotten to know a fair amount about NDSolve from these tutorials:

• NDSolvePlugIns
• NDSolveDesign
• NDSolveStateData

You can also spelunk part of the code for various methods if you root around NDSolve`*:

?NDSolve`*Rung*a

?NDSolve`LSODA

?NDSolve`*inearly*

[etc.]

LSODA is based on LSODE which can be googled. I've found Fortran and C implementations. Obviously, LSODA would have to be rewritten for Mathematica, if for no other reason than to handle arbitrary working precision numbers. Similarly for IDA, which is part of Sundials.

From the tutorials, you find out that most NDSolve methods are implemented in an object-oriented way. The constructor is usually spelunkable and has the form:

LSODA /: NDSolve`InitializeMethod[LSODA,....

or generally

id /: NDSolve`InitializeMethod[id,....

The object created has the form id[data] and the methods have the form id[data][meth] such as id[data]["Step"[..args..]] to advance integration one step.

FEM is a recent extension of NDSolve that is well documented and fairly spelunkable (root around NDSolve`FEM`*), though mesh generation is probably in C. It also relies on LAPACK, which can be googled. See:

• FEM User Guide

I have not yet figured out how to output my own dense output (InterpolationOrder -> All) using the NDSolve framework. (NDSolve will do it for you, using its own algorithm, and I can do it by explicitly constructing an InterpolatingFunction directly. I cannot get NDSolve to use the data I compute to construct the IF, though.)

P.S. You can add TraceInternal -> True to TraceScan in to traceView2[], if you want an organized display of the available internal steps. Beware, you might have to kill Mathematica if the output is of the order of gigabytes. I have had to do so. So do it advisedly. Limit the number of steps NDSolve can take, and don't do it on a huge PDE.