Package for differential equation solvers

AB a at b.c
Sun Jan 10 18:00:04 UTC 2021


On Saturday, 2 January 2021 at 09:47:28 UTC, Flávio wrote:
> Hi folks,
>
> I have started a repository to port Differential equation 
> solver code to D.
> The goal is to have efficient implementations in D that have a 
> much simpler interface for users than current C, C++, and 
> FORTRAN libraries.
>
> Currently is just a couple of solvers quickly whipped together 
> (probably with bugs) as a proof of concept. Please join me to 
> help create this!
>
> https://github.com/fccoelho/D-DifferentialEquations
>
> Contributors are very, very, VERY MUCH welcome!

Hi Flavio, thanks for starting this discussion, it is some time 
that I think D has potential for numeric code, but every time I 
tried to use it I have found that it lacks some features I want.

First, I agree with Kirill that a templated interface is better 
than one based on real slices.
In my experience it is quite useful to
- use higher or lower precision types to trade speed for accuracy 
or vice versa,
- use arbitrary precision types or dynamic precision,
- use auto differentiation on parameters,
- allow matrix(or tensor)-valued equations.

It seems to me that an explicit step can cover the full 
generality by being templated over three  types:
-the type of time
-the type of f(t)
Probably it is possible to deduce all types from the signature of 
f, but that would require a meta-programming expert

Example 1:

auto explicitEulerStep(fun, time, val)( fun f, time t, time dt, 
val y)
// some constraint magic to check that the following operations 
are possible:
// 1) "t+dt"
// 2) "f(t,y)" is a valid call
// 3) "f(t,y)" returns a val
{
return tuple(t+dt, y+dt*f(t,y));
}

Note that there is no need to pass additional parameters to f. 
This is because it can be a delegate (or function literal) and it 
can contains its own parameters.

Example 2: explicitEulerStep( (real t, real y)=>sin(k*t)*y,t, y );

I am a bit less confident about a viable implementation strategy 
for implicit methods.
In this case a solution strategy for the resulting equation is 
needed.
For ease of use a default should be provided, for instance 
findRoot from std.numeric, but that would not work for vector 
valued equations and in general the ability to tailor the solver 
is needed.

auto implicitEulerStep(solver, fun, time, val)( fun f, time t, 
time dt, val y)
{
return tuple(t+dt, solver.find_root( (val yy)=> 
yy-y-dt*f(t+dt,yy)) );
}

Hopefully it is possible to find a nice solution.

Happy new year!

AB



More information about the Digitalmars-d-announce mailing list