non-Newtonian calculus for ODE Solvers

[RobinSchmidt]

after seeing these two videos yesterday:

https://www.youtube.com/watch?v=shdK9DAiDBE
https://www.youtube.com/watch?v=quM45lJUwc4

i digged a little deeper and found some interesting new (to me) math:

https://hal.archives-ouvertes.fr/file/index/docid/945788/filename/nncam.pdf

which could be useful for ODE solvers. i already have some general code for doing stuff like euler-steps, runge-kutta-steps, etc. - but i think they can be straightforwardly generalized like that:

https://www.hindawi.com/journals/aaa/2015/594685/
edit: the link seems to be dead now but i think it was this paper (if i remember correctly):
https://www.researchgate.net/publication/268157925_Generalized_Runge-Kutta_Method_with_respect_to_the_Non-Newtonian_Calculus

...maybe for some problems, these generalized methods yield superior performance over the classical methods? the first paper (which deals mostly with alpha=exp, beta=id (or vice versa, i'm too lazy to check) in the more general context of the 2nd paper) says (page 5, bottom) that the new Euler method based on that will give exact results for the linear differential equation y' = k*y ...that could be a desirable feature

[RobinSchmidt]

i think, the regular Euler-step for the ODE y' = f(y) with stepsize h:
y[n+1] = y[n] + h * f(y[n])
would just translate to:
y[n+1] = y[n] * exp(h * f(y[n]) / y[n])
in what they call the "geometric calculus". this is actually quite similar to what i'm currently doing with the non-uniform filters. they also use an exponential function involving the step-size

[RobinSchmidt]

here's is an interesting article about this starting from the definition of the "multiplicative integral":
https://www.cantorsparadise.com/a-new-kind-of-calculus-5cc5ed84506d

[elanhickler]

Pages and pages of beautiful chaos algorithms! http://www.3d-meier.de/tut19/Seite1.html

Sorry to distract you again, but could you choose any chaos algorithm and make an example using ODE solver? Some of these chaos algorithms break at somewhat disappointingly low frequencies (due to step size of higher frequencies). Maybe some better interpolation will make them run better.

[RobinSchmidt]

Oh wow! Thanks! I've just seen your post today. That's indeed a nice collection! So - you have implemented a simple numerical solver like forward Euler and are interested in trying out better ODE solver algorithms on these? I'm not sure, if this non-Newtonian approach will help here but I also have some experimental code for higher order ODE solving methods somewhere in the prototypes and experiments... Runge-Kutta, Adams-Bashforth, Adams-Moulton, etc.

[elanhickler]

forward Euler? No, I've just been using the naive equations. I'll look for the higher order ODE solves in the library.

[RobinSchmidt]

Ah - damn! There is not yet much in there yet. I once started writing this stub:

https://github.com/RobinSchmidt/RS-MET/blob/master/Libraries/RobsJuceModules/rs_testing/Prototypes/OdeSolver.h

and I have some old, unused legacy code here:

https://github.com/RobinSchmidt/RS-MET/blob/master/Misc/UnusedCode/RSLibAndTests/RSLib/Code/Math/Misc/DifferentialEquationSystem.h

on which I want to (partially) base it. .....but only formula/algorithm-wise, not API-wise because the old legacy API was horribly inefficient (it was meant to be a prototype only anyway). I actually like to have efficient explicit Runge-Kutta and Adams-Bashforth solvers at some point. And maybe later also implicit solvers.

Edit: In the GNUPlotCpp codebase, there's also a simple solver:

https://github.com/RobinSchmidt/GNUPlotCPP/blob/master/Code/MathTools.h

This is one with stepsize control - probably not so suitable for a realtime signal generator, but needed for creating some plots. I've used it to produce these plots:

RobinSchmidt/GNUPlotCPP#1

[RobinSchmidt]

I think, for a realtime application, it may be a good idea to compute an initial section using Runge-Kutta but when enough samples have been computed, switch over to an Adams-Bashforth method. The latter is a multistep method which needs several past samples to compute the next - and for the initial section, you don't have any so you have to compute them with a different method like Runge-Kutta. That would probably give a nice efficient realtime solver...or so I hope

[RobinSchmidt]

If you are just doing something like:

yp = f(y);     // compute derivative: y' = f(y)
y += h*yp;     // do update step: yNew = yOld + h * y'

then that is forward Euler. The simplemost ODE solver technique in existence.

[elanhickler]

class Attractor
{
    void setSampleRate(double newSampleRate)
    {
        sampleRate = newSampleRate;
        update();
    }

    void setFrequency(double v)
    {
        frequency = v;
        update();
    }

    /** Sets the internal state of the system consisting of the 3 state variables x, y, z. */
    void setState(double inX, double inY, double inZ)
    {
        x = inX;
        y = inY;
        z = inZ;
    }

    /* Because each attractor has a different set of constants, we use an array to set values arbitrarily */
    void setConstant(uint index, double value)
    {
        C[index] = value;
    }

    double getX()
    {
       return x;
    }

    double getY()
    {
       return y;
    }

    double getZ()
    {
       return z;
    }

    void inc(){}

    void reset(){}

    protected void update()
    {
        h = frequency / sampleRate;
    }    

    protected array<double> C(10); // array of constants
    protected double sampleRate = 44100;
    protected double frequency = 100;
    protected double h = frequency / sampleRate;
    protected double x, y, z, w;
    protected double dx, dy, dz, dw;
}

class Lorentz : Attractor
{
    uint sigma = 0;
    uint rho   = 1;
    uint beta  = 2;

    Lorentz()
    {        
        C[sigma] = 10.0;
        C[rho]   = 28.0;
        C[beta]  = 8.0/3.0;

        reset();
    }

    void inc()
    {
        dx = C[sigma] * (y - x);
        dy = x * (C[rho] - z) - y;
        dz = x * y - C[beta] * z;

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = 0.5;
        y = 0.0;
        z = 0.0;
    }
};

class Aizawa : Attractor // http://www.3d-meier.de/tut19/Seite3.html
{
    uint a = 0;
    uint b = 1;
    uint c = 2;
    uint d = 3;
    uint e = 4;
    uint f = 5;

    Aizawa()
    {
        C[a] = 0.95;
        C[b] = 0.7;
        C[c] = 0.6;
        C[d] = 3.5;
        C[e] = 0.25;
        C[f] = 0.1;

        reset();
    }

    void inc()
    {
        dx = (z - C[b]) * x - C[d] * y;
        dy = C[d] * x + (z - C[b]) * y;
        dz = C[c] + C[a] * z - (pow(z,3) / 3.0) - (x*x + y*y) * (1 + C[e] * z) + C[f] * z * pow(x,3);

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = .1;
        y = 0;
        z = 0;
    }
}

class Bouali : Attractor // http://www.3d-meier.de/tut19/Seite5.html
{
    uint a = 0;
    uint s = 1;

    Bouali()
    {
        C[a] = 0.3;    
        C[s] = 1.0;

        reset();
    }

    void inc()
    {
        dx = x * (4-y) + C[a] * z;
        dy = -y * (1 - x * x);
        dz = -x * (1.5 - C[s] * z) -0.05 * z;

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = 1;
        y = 0.1;
        z = 0.1;
    }
};

class ChenLee : Attractor // http://www.3d-meier.de/tut19/Seite8.html
{
    uint a = 0;
    uint b = 1;
    uint c = 2;

    ChenLee()
    {
        C[a] = 5;
        C[b] = -10;
        C[c] = -0.38;

        reset();
    }

    void inc()
    {
        dx = C[a] * x - y * z;
        dy = C[b] * y + x * z;
        dz = C[c] * z + x*y/3.0;

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = 1;
        y = 0;
        z = 4.5;
    }
}	

class DequanLi : Attractor // http://www.3d-meier.de/tut19/Seite9.html
{
    uint a = 0;
    uint c = 1;
    uint d = 2;
    uint e = 3;
    uint k = 4;
    uint f = 5;

    DequanLi()
    {
        C[a] = 40;
        C[c] = 1.833;
        C[d] = 0.16;
        C[e] = 0.65;
        C[k] = 55;
        C[f] = 20;

        reset();
    }

    void inc()
    {
        dx = C[a] * (y - x) + C[d] * x * z;
        double dy = C[k] * x + C[f] * y - x * z;
        dz = C[c] * z + x * y - C[e] * pow(x,2);

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = .349;
        y = 0;
        z = -.16;
    }
}

class DenTSUCS2 : Attractor // http://www.3d-meier.de/tut19/Seite43.html
{
    uint a = 0;
	uint b = 1;
    uint c = 2;
    uint d = 3;
    uint e = 4;
    uint f = 5;

    DenTSUCS2()
    {
        C[a] = 40;
		C[b] = 55;
        C[c] = 1.833;
        C[d] = 0.16;
        C[e] = 0.65;
        C[f] = 20;

        reset();
    }

    void inc()
    {
        dx = C[a] * (y - x) + C[d] * x * z;
        dy = C[b] * x - x * z + C[f] * y;
        dz = C[c] * z + x * y - C[e] * pow(x,2);

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = .1;
        y = 1;
        z = -.1;
    }
}

class DenGenesioTesi : Attractor // http://www.3d-meier.de/tut19/Seite11.html
{
    uint a = 0;
	uint b = 1;
    uint c = 2;

    DenGenesioTesi()
    {
        C[a] = .44;
		C[b] = 1.1;
        C[c] = 1.0;

        reset();
    }

    void inc()
    {
        dx = y;
        dy = z;
        dz = -C[c] * x - C[b] * y - C[a] * z + pow(x,2);

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = .1;
        y = .1;
        z = .1;
    }
}

class Hadley : Attractor // http://www.3d-meier.de/tut19/Seite12.html
{
    uint a = 0;
	uint b = 1;
	uint f = 2;
    uint g = 3;

    Hadley()
    {
        C[a] = .2;
		C[b] = 4;
        C[f] = 8;
		C[g] = 1;

        reset();
    }

    void inc()
    {
        dx = -pow(y,2) - pow(z,2) - C[a] * x + C[a] * f;
        dy = x*y - C[b]*x*z - y + C[g];
        dz = C[b]*x*y + x*z - z;

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = .1;
        y = 0;
        z = 0;
    }
}

class Halvorsen : Attractor // http://www.3d-meier.de/tut19/Seite13.html
{
    uint a = 0;

    Halvorsen()
    {
        C[a] = 1.4;

        reset();
    }

    void inc()
    {
        dx = -C[a]*x - 4*y - 4*z - pow(y,2);
        dy = -C[a]*y - 4*z - 4*x - pow(z,2);
        dz = -C[a]*z - 4*x - 4*y - pow(x,2);

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = 1;
        y = 0;
        z = 0;
    }
}

class HyperchaoticQi : Attractor // http://www.3d-meier.de/tut19/Seite90.html
{
    uint a = 0;	
    uint b = 1;
    uint c = 2;
    uint d = 3;
    uint e = 4;
    uint f = 5;

    HyperchaoticQi()
    {
        C[a] = 50;
        C[b] = 24;
        C[c] = 13;
        C[d] = 8;
        C[e] = 33;
        C[f] = 30;

        reset();
    }

    void inc()
    {
        dx = C[a] * (y - x) + y * z;
        dy = C[b] * (x + y) - x * z;
        dz = -C[c]*z - C[e]*w + x*y;
        dw = -C[d]*w + C[f]*z + x*y;

        x += h * dx;
        y += h * dy;
        z += h * dz;
        w += h * dw;
    }

    void reset()
    {
        x = 10;
        y = 15;
        z = 20;
        w = 22;
    }
}

class Dadra : Attractor //http://www.3d-meier.de/tut19/Seite77.html
{
    uint p = 0;	
    uint q = 1;
    uint r = 2;
    uint s = 3;
    uint e = 4;

    Dadra()
    {
        C[p] = 3;
        C[q] = 2.7;
        C[r] = 1.7;
        C[s] = 2;
        C[e] = 9;

        reset();
    }

    void inc()
    {
        dx = y - C[p]*x + C[q]*y*z;
        dy = C[r]*y - x*z + z;
        dz = C[s]*x*y - C[e]*z;

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = .1;
        y = .1;
        z = .1;
    }
}

class DenSprottLinzJ : Attractor // http://www.3d-meier.de/tut19/Seite29.html
{
    uint a = 0;	

    DenSprottLinzJ()
    {
        C[a] = 2;

        reset();
    }

    void inc()
    {
        dx = C[a]*z;
        dy = -C[a]*y + z;
        dz = -x + y + y*y;

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = .1;
        y = .1;
        z = .1;
    }
}

class LiuChen : Attractor // http://www.3d-meier.de/tut19/Seite46.html
{
    uint a = 0;
    uint b = 1;
    uint c = 2;
    uint d = 3;
    uint e = 4;
    uint f = 5;
    uint g = 6;    

    LiuChen()
    {
        C[a] = 2.4;
        C[b] = -3.78;
        C[c] = 14;
        C[d] = -11;
        C[e] = 4;
        C[f] = 5.58;
        C[g] = -1;

        reset();
    }

    void inc()
    {
        dx = C[a]*y + C[b]*x + C[c]*y*z;
        dy = C[d]*y - z + C[e]*x*z;
        dz = C[f]*z + C[g]*x*y;

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = 1;
        y = 3;
        z = 5;
    }
}

class MultiSprottC : Attractor // http://www.3d-meier.de/tut19/Seite192.html
{
    uint a = 0;
    uint b = 1;

    MultiSprottC()
    {
        C[a] = 2.5;
        C[b] = 1.5;

        reset();
    }

    void inc()
    {
        dx = C[a] * (y - x);
        dy = x * z;
        dz = C[b] - y*y;

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = .1;
        y = .1;
        z = .1;
    }
}

class SprottLinz : Attractor // http://www.3d-meier.de/tut19/Seite20.html
{
    SprottLinz()
    {
        reset();
    }

    void inc()
    {
        dx = y;
        dy = -x + y*z;
        dz = 1 - y*y;

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = .1;
        y = .1;
        z = .1;
    }
}

class Thomas : Attractor // http://www.3d-meier.de/tut19/Seite20.html
{
    uint b = 0;

    Thomas()
    {
        C[b] = .19;

        reset();
    }

    void inc()
    {
        dx = -C[b]*x + sin(y);
        dy = -C[b]*y + sin(z);
        dz = -C[b]*z + sin(x);

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = .1;
        y = 0;
        z = 0;
    }
}

class LorenzMod1  : Attractor // http://www.3d-meier.de/tut19/Seite79.html
{
    uint a = 0;
    uint b = 1;
    uint c = 2;
    uint d = 3;

    LorenzMod1()
    {
        C[a] = 0.1;
        C[b] = 4;
        C[c] = 14;
        C[d] = 0.08;

        reset();
    }

    void inc()
    {
        dx = -C[a]*x + y*y - z*z + C[a]*C[c];
        dy = x * (y - C[b]*z) + C[d];
        dz = z + x * (C[b]*y + z);

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = .1;
        y = .1;
        z = .1;
    }
}

class LorXYZ15  : Attractor // https://arxiv.org/ftp/arxiv/papers/1409/1409.7842.pdf
{
    uint a = 0;
    uint b = 1;
    uint c = 2;

    LorXYZ15()
    {
        C[a] = 2;
        C[b] = 0.82942;
        C[c] = 7.178;

        reset();
    }

    void inc()
    {  
        double sinY = sin(y);

        dx = z + sin(pow(y,C[a])) * (z + sinY);
        dy = z + sin(x*y) * (sinY - C[b]);
        dz = x * (sinY - C[c]);

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = .1;
        y = .1;
        z = .1;
    }
}

class WangSun  : Attractor // http://www.3d-meier.de/tut19/Seite99.html
{
    uint a = 0;
    uint b = 1;
    uint c = 2;
    uint d = 3;
    uint e = 4;
    uint f = 5;

    WangSun()
    {
        C[a] = .2;
        C[b] = -.01;
        C[c] = 1;
        C[d] = -0.4;
        C[e] = -1;
        C[f] = -1;

        reset();
    }

    void inc()
    {  
        dx = x*C[a] + C[c]*y*z;
        dy = C[b]*x + C[d]*y - x*z;
        dz = C[e]*z + C[f]*x*y;

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = .3;
        y = .1;
        z = 1;
    }
}

class HindenmarshRose : Attractor // https://en.wikipedia.org/wiki/Hindmarsh%E2%80%93Rose_model
{
    uint a = 0;
    uint b = 1;
    uint c = 2;
    uint d = 3;
    uint r = 4;
    uint xR = 5;
    uint s = 6;
    uint I = 7;

    HindenmarshRose()
    {
        C[a] = 1;
        C[b] = 3;
        C[c] = 1;
        C[d] = 5;
        C[r] = .005;
        C[s] = 4.0;
        C[xR] = -1.6;
        C[I] = 3.25;

        reset();
    }

    void inc()
    {  
        dx = y + (C[b] * x * x) - (C[a] * x * x * x) - z + C[I];
        dy = C[c] - (5 * x * x) - y;
        dz = C[r] * (C[s] * (x - C[xR]) - z);

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

    void reset()
    {
        x = -1.6; // membrane potential
        y = 4.0;
        z = 2.75;
    }
}

[RobinSchmidt]

Nice!

yep - that is precisely forward Euler:

    void inc()
    {
        dx = C[sigma] * (y - x);
        dy = x * (C[rho] - z) - y;
        dz = x * y - C[beta] * z;

        x += h * dx;
        y += h * dy;
        z += h * dz;
    }

You just compute the derivatives x', y', z' via the given formulas and then take a step with stepsize h into that computed direction

[RobinSchmidt]

OK - so which one of these is a good example for breaking at disappointingly low frequencies? Is it OK, if I copy this code into my research codebase for some experimentation?

[elanhickler]

yes, this is free to use. I was planning on translating a lot more chaos algorithms and making a VST plugin, a "chaos algorithm explorer" or something.

class DenTSUCS2 : Attractor // http://www.3d-meier.de/tut19/Seite43.html

Note: This one runs loud.

[RobinSchmidt]

I was planning on translating a lot more chaos algorithms

After adding the "Formula" module to Liberty, I had the idea of making a similar module for an ODE solver. Something that just let's you enter the 3 formulas (or 2 or 4 or whatever number the system has) and then runs an ODE solver using these formulas in an expression evaluator.

[RobinSchmidt]

OK - nice. This is fun. Some first results are in:

The code is here - in function testAttractors, currently (2022/06/19) starting at line 9104 (it's a big messy file with all sorts of experiments):

https://github.com/RobinSchmidt/RS-MET-Research/blob/master/Prototypes/Cpp/Source/Experiments.cpp

I have used the forward Euler method that was already there and I also added a function that performs steps according to the midpoint rule. This did indeed increase the usable frequency range by almost an order of magnitude (factor 9 or something). Perhaps even higher order methods would extend that even more. ...more tests needed....but it's probably a good idea to implement these methods in a generic ODE solver class rather than coding the same math ad hoc over and over again for every system separately. That kind of code duplication may be ok for a very simple method like forward Euler and maybe midpoint, too. Maybe like borderline acceptable for the latter - but already for something just a little bit more complex (like, say, 4th order Runge-Kutta - which is really just the beginning when it comes to ODE methods), it's really nicer to get rid of the code duplication and implement the method just generically once and for all and then call it for the particular system at hand through a sensible API. I think, this may be a perfect use case for std::function.

btw: don't worry about being "loud". I'm just staring at plots. :-)

[RobinSchmidt]

Oh and your attractor code is here:

https://github.com/RobinSchmidt/RS-MET-Research/blob/master/Prototypes/Cpp/Source/Attractors.cpp

I needed to tweak it quite a bit to make it compile. There were semicolons missing, public declarations for the inheritance, etc. What compiler are you using that accepted the original code? Is this allowed in C++17 or 20 or something?

Oh - by the way: I actually think, solvers based this non-Newtonian calculus idea would be suitable only for strictly positive functions, since there appears this division by y in my (proposed) formula for an adapted Euler rule....but maybe I'm on the wrong track with my proposed formula anyway...don't really know for sure yet