Why does using a smaller tolerance for error result in a worse estimate when using Runge-Kutta Fehlberg method?

Question

I developped an ODE solver which uses Fehlberg's Adaptive Step Size procedure, but it appears that when I go to very small tolerances that the "Optimal" step size becomes very large and then my final estimate is very far from the true value of the function (worse than regular runge-kutta)

I was wondering if anyone would know why this is and how I could remedy it. I am using the Cash-Karp parameters and the following equations to calculate new step size ('h' is the step size):

    delta0 = error * (Abs(y4old) + Abs(h * MyODE(x, y4old)))   
    delta1 = Abs(y5 - y4)

    If delta1 > error Then
        If delta0 > delta1 Then
            h = 0.9 * h * (Abs(delta1 / delta0)) ^ 0.2   '0.9 is a safety factor
        Else
            h = 0.9 * h * (Abs(delta1 / delta0)) ^ 0.25
        End If
    End If

Here is the full code and an example of what I mean:

Public Function Fehlberg(x As Double, xmax As Double, y4 As Double, h As Double, error As Double) As Double

Dim k(7) As Double, cond As Boolean, delta0 As Double, delta1 As Double, y4old As Double

While x < xmax

    If x + h < xmax Then
        x = x + h
    Else
        h = xmax - x
        x = xmax
    End If

    k(1) = MyODE(x, y4)
    k(2) = MyODE(x + h * 0.2, y4 + h * k(1) * 0.2)
    k(3) = MyODE(x + 0.3 * h, y4 + h * (0.075 * k(1) + 0.225 * k(2)))
    k(4) = MyODE(x + 0.6 * h, y4 + h * (0.3 * k(1) - 0.9 * k(2) + 1.2 * k(3)))
    k(5) = MyODE(x + h, y4 + h * ((-11 / 54) * k(1) + 2.5 * k(2) - (70 / 27) * k(3) + (35 / 27) * k(4)))
    k(6) = MyODE(x + 0.875 * h, y4 + h * ((1631 / 55296) * k(1) + (175 / 512) * k(2) + (575 / 13824) * k(3) + (44275 / (110592) * k(4) + (253 / 4096) * k(5))))

    y4old = y4                                          'Saving old y4 value to calc delta0
    y4 = y4 + h * (k(1) * (37 / 378) + k(3) * (250 / 621) + k(4) * (125 / 594) + k(6) * (512 / 1771))
    y5 = y4 + h * (k(1) * (2825 / 27648) + k(3) * (18575 / 48384) + k(4) * (13525 / 55296) + k(5) * (277 / 14336) + k(6) * (0.25))

    delta0 = error * (Abs(y4old) + Abs(h * MyODE(x, y4old)))    'or just error * y
    delta1 = Abs(y5 - y4)

    If delta1 > error Then
        If delta0 > delta1 Then
            h = 0.9 * h * (Abs(delta1 / delta0)) ^ 0.2   '0.9 is a safety factor
        Else
            h = 0.9 * h * (Abs(delta1 / delta0)) ^ 0.25
        End If
    End If
Wend

Fehlberg = y4

End Function

So if my ODE is y'=2x (y = x^2 + C) and I want to find y(3) starting at y(0) = 1, with a tolerance of 0.01 I get an answer of 10.02544 but with a tolerance of 0.001 (which should give closer to the true value) I get 17.76559251. Why is this happening?

Thank you in advance for any help you provide!

This may be the wrong place to post this question so please let me know where I should post if not here.

Answer 1

As you increase x before computing the steps, you have to modify your step computation to use x-h as base point, i.e.,

k(1) = MyODE(x-h, y4)
k(2) = MyODE(x-h + h * 0.2, y4 + h * k(1) * 0.2)

or

k(2) = MyODE(x - h * 0.8, y4 + h * k(1) * 0.2)

to apply the method according to its theory. Doing that returns the correct result for the test example independent of the error threshold, as any order 4 method should for polynomials in x up to degree 4.

Applying the same for the example $y'=y$ leads again to strange behavior indicating that the step size selector is not well thought-out. One error is that the fraction is the wrong way, as delta0 is the desired error and delta1 is the estimated error, the formulas should be

h = 0.9 * h * (delta0/delta1)^0.25

etc.

---

A more correct stepping method should lead to behavior like in the table

input      approx. of exp(3)   difference
error     20.0855369231876679
-------------------------------------------
1.00e-02  20.010411826192708  -7.5125e-02 
1.00e-03  20.010411826192708  -7.5125e-02 
1.00e-04  20.079949162409612  -5.5878e-03 
1.00e-05  20.085081008767098  -4.5591e-04 
1.00e-06  20.085497474964928  -3.9448e-05 
1.00e-07  20.085533801160185  -3.1220e-06 
1.00e-08  20.085536648684410  -2.7450e-07 
1.00e-09  20.085536898129003  -2.5059e-08 
1.00e-10  20.085536920778321  -2.4093e-09 
1.00e-11  20.085536922954645  -2.3302e-10 
1.00e-12  20.085536923164945  -2.2723e-11 
1.00e-13  20.085536923185341  -2.3270e-12 
1.00e-14  20.085536923187494  -1.7408e-13 
1.00e-15  20.085536923187721   5.3291e-14 

so that indeed the relative error of the result is comparable to the given tolerance.