2
$\begingroup$

I have been working on creating a few home-made numerical methods, and I am using them to visualize text-book problems from my Strogatz dynamics textbook. It feels like a good way to learn numerical methods and dynamics at the same time.

A recent problem from the book deals with "Critical Slowdown" and is exemplified by comparing $\dot{x} = -x^3$ with $\dot{x} = -x$. I can't pretend to understand anything except that one decays more quickly than the other.

I decided to run both through my little home-made RK4 algorithm and graph them. What I found is that I get overflow for too-large step size, and I was wondering why. Conceptually, what causes this type of behavior?

I know that the step size should be chosen in accordance with how quickly the solution varies, or error will be very high. What I don't get is the overflow errors I'm getting.

Code (The ODE's are defined and solved at the bottom. To see the overflow, replace .001 with .1 in the step-size argument.:

def rk4(dt, t, field, y_0):

    y = np.asarray(len(t) * [y_0]) 

    for i in np.arange(len(t) - 1):
        k1 = dt * field(t[i], y[i])
        k2 = dt * field(t[i] + 0.5 * dt, y[i] + 0.5 * k1)
        k3 = dt * field(t[i] + 0.5 * dt, y[i] + 0.5 * k2)
        k4 = dt * field(t[i] + dt, y[i] + k3)
        y[i + 1] = y[i] + (k1 + 2 * k2 + 2 * k3 + k4) / 6

    return y

def vf_grapher(fn, t_0, t_n, dt, y_0, lintype='-r', sup_title=None,
               title=None, xlab=None, ylab=None):

    t = np.arange(t_0, t_n, dt)
    y_min = .0
    y_max = .0

    fig, axs = plt.subplots()
    fig.suptitle(sup_title)

    axs.set_title(title)
    axs.set_ylabel(ylab)
    axs.set_xlabel(xlab)

    for iv in np.array(y_0, ndmin=1, copy=False):
        soln = rk4(dt, t, fn, iv)
        plt.plot(t, soln, lintype)
        if y_min > np.min(soln):
            y_min = np.min(soln)
        if y_max < np.max(soln):
            y_max = np.max(soln)

    x = np.linspace(t_0, t_n + dt, 11)
    y = np.linspace(y_min, y_max, 11)

    X, Y = np.meshgrid(x, y)

    theta = np.arctan(fn(X, Y))

    U = np.cos(theta)
    V = np.sin(theta)

    plt.quiver(X, Y, U, V, angles='xy')
    plt.xlim((t_0, t_n - dt))
    plt.ylim((y_min - .05 * y_min, y_max + .05 * y_max))
    plt.show()


if __name__ == '__main__':

    def crit_slow(t, x):
        return -x**3


    def exp_decay(t, y):
        return -y

    # graph an ODE modelling critical slow down (has discontinuities in it's 
    second derivative)
    vf_grapher(crit_slow, 0., 10., .001, 10.)

    # graph an ODE modelling exponential decay (compare to critical slowdown)
    vf.vf_grapher(exp_decay, 0., 10., .001, 10.)
$\endgroup$
  • 1
    $\begingroup$ All explicit time integration schemes have conditional stability, meaning they will avoid blowing up (and overflowing) as long as you keep your timestep below some upper bound. In the case of nonlinear problems, this upper bound can change from timestep to timestep, so that brings in extra complexity. So in fact, choosing a scheme to use can depend on your stability, accuracy, and potentially even computational complexity requirements. $\endgroup$ – spektr Oct 31 '19 at 16:25
  • $\begingroup$ @spektr does that stability have to do with so-called CFL requirements? Or am I getting my reading mixed up? $\endgroup$ – rocksNwaves Oct 31 '19 at 16:41
  • $\begingroup$ That’s on the right track, except the CFL condition is primarily for PDEs. Each explicit time integration scheme (eg explicit Euler, RK4, etc) inherently require a bound on their stepsize to ensure the errors don’t “blow up”, which generally results in overflow and all that terrible stuff. If you look through a numerical analysis book on these methods, they should at least tell you the stability region, if not also prove the bound. $\endgroup$ – spektr Oct 31 '19 at 16:44
  • $\begingroup$ @spektr, Thank you. I just ordered a few textbooks on ODEs and Dynamics. Can you recommend one for Numerical Analysis? In particular, understandable by one who is slow on the uptake. I just graduated with a master's in applied math but I feel like I didn't learn anything and was given a passing grade out of sympathy. Now I'm trying to ACTUALLY learn this stuff. $\endgroup$ – rocksNwaves Oct 31 '19 at 16:53
  • 1
    $\begingroup$ One that’s pretty easy to follow (though you may want something more as an applied math grad) is Fundamentals of Engineering Numerical Analysis, Second Edition by Moin. It goes over more than just solving differential equations, but it does cover the material on stability. Consider trying to find a PDF of the book before buying it because you may only need a small part of the text for your purposes. $\endgroup$ – spektr Oct 31 '19 at 16:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.