Problem with Heun's solver in 'desolver'

[tp-Euler]

Dear Ekin,

I appreciate your Python package 'desolver'. In this package, the methods 'Euler', 'RK4' or 'Runge-Kutta 4' seem to work fine.

However, I am writing you to inform you that the Heun method (see e.g. https://en.wikipedia.org/wiki/Heun's_method) does not seem to be properly implemented. For that, I did try all of the method's which could potentially represent Heun's method, i.e.

'Euler-Trap', 'Euler-Trapezoidal', 'EulerSolver', 'EulerTrapSolver', 'Explicit Euler-Trapezoidal', "Explicit Heun's", 'Explicit Midpoint', "Heun's", 'HeunEulerSolver', 'HeunsSolver',

And which of these method is in fact supposed to be the Heun method, i.e. the one described at the above wikipedia link?

More precisely, I did use your package 'desolver' in order to solve the ODE dx/dt = - sqrt(x/t) with initial condition x(t=1.0) = x_0 = 1.0 and stepsize of Dt=0.1. With Heun's method (as described on the wikipedia link above), one obtains
x_1 = x(t=1.1) = 0.9047733 and
x_2 = x(t=1.2) = 0.81824429
But I do not get these values with any of the methods listed above. By the way, some of the methods listed above also yield the very same values (for x_1 and x_2).

Could you please look into this soon?

Thank you and best regards,
Thomas

[Microno95]

Hi Thomas,

That's very odd and interesting! All the methods do pass the unit tests, so maybe I need to include a few more equations to test.

Could you provide a small snippet of code that reproduces the issue? In addition, please let me know the exact version of the desolver library you're using, the platform you're running on and the versions of Python, numpy, pytorch (if you're using the torch backend) and jax (if you're using the jax backend). I will try and reproduce the issue and get back to you with my results.

The Heun's method described in wikipedia is all of ["Explicit Heun's", "Heun's", "HeunsSolver"], they are basically aliases to the same solver underneath. Heun-Euler is a separate adaptive method that uses Heun's for the 2nd order and Euler for the 1st order estimate as part of an embedded Runge-Kutta scheme.

Best,
Ekin

[Microno95]

I'll also note that "Heun's Method" is a name given to two Runge-Kutta methods: Improved Euler and Ralston's Method.

[tp-Euler]

Hi Ekin,

I use version 5.0.0 of desolver.

And here's a snippet of code that reproduces the issue:

import numpy as np
import desolver as de

# RHS: x' = -sqrt(x/t) 
def rhs(t, x, **kwargs):
    return -np.sqrt(x/t)

x_0 = np.array([1.0])

dgl = de.OdeSystem(rhs, y0=x_0, t=(1., 2.), dt=0.1)

# aktuelle Methode (Default)
print("default method:", dgl.method)

# verfügbare Methoden anzeigen (gute Möglichkeit, den richtigen Namen zu finden)
print("available methods:", de.integrators.available_methods(names=True))

# Methode auf 'Euler-Cauchy' setzen — Name kann je nach Version auch 'heun' o.ä. sein
dgl.set_method('HeunsSolver')   

dgl.integrate()

# extract results
t = dgl.t   # array bestehend aus den Zeitpunkten
x = dgl.y   # für 1D-Problem: y ist das Array der y-Werte

With the code above, I get
x_1 = 0.90484392 and
x_2 = 0.8183645
rather than the values in my first message above.

Best, Thomas

[Microno95]

(I have wrapped the code snippet in "```" which treats it as a code block)

I will test this and get back to you.

[Microno95]

Hi Thomas,

Having tested your code, I have concluded that what you want is the Heun's Method that is not Ralston's Method (the one implemented in DESolver).

The following code snippet gets you the exact integrator you want and numerically outputs the exact values you expect:

import numpy as np
import desolver as de

x_0 = np.array([1.0])

dgl = de.OdeSystem(rhs, y0=x_0, t=(1., 2.), dt=0.1)

class HeunsFromWikipedia(de.integrators.RungeKuttaIntegrator):
    """
    The derived class that implements Heun's method.
    """
    
    __order__ = 2.0

    __alt_names__ = ("Explicit Heun's", "Heun's")

    tableau_intermediate = np.array(
        [[0,    0,   0  ],
         [1,    1,   0  ]], dtype=np.float64
    )

    tableau_final = np.array(
        [[0,    1/2, 1/2]], dtype=np.float64
    )

# Methode auf 'Euler-Cauchy' setzen — Name kann je nach Version auch 'heun' o.ä. sein
dgl.set_method(HeunsFromWikipedia)

dgl.integrate()

# extract results
t = dgl.t   # array bestehend aus den Zeitpunkten
x = dgl.y   # für 1D-Problem: y ist das Array der y-Werte

for idx, (current_time, current_state) in enumerate(zip(t, x)):
    print(f"x(t={current_time:.4f}) {current_state[0]:.16f}")

Outputs:

x(t=1.0000) 1.0000000000000000
x(t=1.1000) 0.9047732983133354
x(t=1.2000) 0.8182442923083457
x(t=1.3000) 0.7393434415151319
x(t=1.4000) 0.6672041133207696
x(t=1.5000) 0.6011124220267992
x(t=1.6000) 0.5404719771992487
x(t=1.7000) 0.4847784552940045
x(t=1.8000) 0.4336008410356030
x(t=1.9000) 0.3865673164150622
x(t=2.0000) 0.3433544623760296

If you do want to sample specific times (ie. t=1.1, t=1.2, etc.), I would recommend that you use event detection and an adaptive integrator which ensures the numerical accuracy of your results. A small change to your code would be to run:

import numpy as np
import desolver as de

# RHS: x' = -sqrt(x/t) 
def rhs(t, x, **kwargs):
    return -np.sqrt(x/t)

x_0 = np.array([1.0])

dgl = de.OdeSystem(rhs, y0=x_0, t=(1., 2.), dt=0.001)

# aktuelle Methode (Default)
print("default method:", dgl.method)

# verfügbare Methoden anzeigen (gute Möglichkeit, den richtigen Namen zu finden)
print("available methods:", de.integrators.available_methods(names=True))

# Methode auf 'Euler-Cauchy' setzen — Name kann je nach Version auch 'heun' o.ä. sein
dgl.set_method('RK87')

def sample_time_event(t, y):
    return 10*t - np.round(10*t)

def constrain_step_size_cb(ode_sys):
    ode_sys.dt = np.minimum(ode_sys.dt, 0.05)

dgl.integrate(callback=constrain_step_size_cb, events=sample_time_event)

# extract results
t = np.stack([i.t for i in dgl.events if i.event == sample_time_event])   # array bestehend aus den Zeitpunkten
x = np.stack([i.y for i in dgl.events if i.event == sample_time_event])   # für 1D-Problem: y ist das Array der y-Werte

for idx, (current_time, current_state) in enumerate(zip(t, x)):
    print(f"x(t={current_time:.4f}) {current_state[0]:.16f}")

Outputs:

x(t=1.0000) 1.0000000000000000
x(t=1.1000) 0.9047645885444331
x(t=1.2000) 0.8182195260720772
x(t=1.3000) 0.7392982890835372
x(t=1.4000) 0.6671361653864403
x(t=1.5000) 0.6010205080328347
x(t=1.6000) 0.5403557386156153
x(t=1.7000) 0.4846380716948372
x(t=1.8000) 0.4334368506042893
x(t=1.9000) 0.3863804963499921
x(t=2.0000) 0.3431457505076204

The key changes here are:

1. We use the sample_time_event function to detect when the time passes a specific value (in this case each tenth of a second),
2. We use the constrain_step_size_cb to ensure that the step size stays within half the sampling time (otherwise we would lose the ability to detect the events) - if you're familiar with sampling theory, this is set by the Nyquist frequency essentially.

Let me know if you have any questions!

[tp-Euler]

Hi Ekin,

I will first stick to your first version as I use it for teaching purposes. Could you biefly tell me, what the line

__alt_names__ = ("Explicit Heun's", "Heun's")

in your above class definition actually does?

I thought these would be alternative names, so that I could replace

dgl.set_method(HeunsFromWikipedia)

e.g. by

dgl.set_method("Explicit Heun's")

but then I again get the previous values. So the method "Explicit Heun's" is actually not overwritten by running the above class definition?

[Microno95]

Hi Thomas,

Ah yes, I see the confusion. I had copied the code from the DESolver internals, but those names are not registered with the internal mapping that maintains the relationship between strings and their actual integrator counterparts which is done on initialisation here:

desolver/desolver/integrators/__init__.py

Lines 50 to 56 in 34602f0

	__available_methods.update(dict(
	[(func.__name__, func) for func in __explicit_integration_methods__ if hasattr(func, "__alt_names__")] +
	[(alt_name, func) for func in __explicit_integration_methods__ if hasattr(func, "__alt_names__") for alt_name in func.__alt_names__]))
	__available_methods.update(dict(
	[(func.__name__, func) for func in __implicit_integration_methods__ if hasattr(func, "__alt_names__")] +
	[(alt_name, func) for func in __implicit_integration_methods__ if hasattr(func, "__alt_names__") for alt_name in func.__alt_names__]))

.

Suffice to say that those lines do nothing in this case. If you'd like to submit a Pull Request with an update to the naming/butcher tableau of the methods in DESolver, I would be more than happy to review that quickly and merge it. Alternatively, while I am currently quite busy, I might be able to look at this in the coming 3 weeks if that suits your timeline.

[tp-Euler]

Hi Ekin,

I actually had just signed up to GitHub yesterday (before writing to you by posting the above "issue") - so I am not familiar with all the collaboration processes on GitHub (such as working with forks, branches, pull requests, etc.)

However, I would appreciate if you could update that naming as discussed above at your earliest convenience.

Many thanks, Thomas

[Microno95]

Hi Thomas, I absolutely understand. I will let you know when that's done. - Ekin