|Posted by AB|
in reply to Flávio
Posted in reply to Flávio
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!
> 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
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!