12
$\begingroup$

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).

Improve this question
$\endgroup$
3
  • $\begingroup$ The article states that the source code of functions written in C cannot be accessed. My guess is that much of NDSolve is in C for efficiency, or that it calls functions from other libraries. Very interesting question, though. $\endgroup$ Commented Mar 13, 2022 at 15:59
  • 2
    $\begingroup$ @bbgodfrey I read somewhere sometime ago that about 50% of Mathematica is written in C and C++ and the rest is in Mathematica (i.e. Wolfram Language). $\endgroup$ Commented Mar 13, 2022 at 16:29
  • $\begingroup$ @bbgodfrey I agree that lower-level NDSolve calls are surely in C/C++, but I'm also confident that there's more of a Wolfram Language backend to NDSolve than what appears here. I ran into the same situation with NIntegrate. $\endgroup$ Commented Mar 13, 2022 at 16:52

1 Answer 1

Reset to default
13
$\begingroup$

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

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:

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Great answer, thank you. $\endgroup$ Commented Mar 13, 2022 at 23:12