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

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?

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

• Could you describe in some greater length the theory or heuristics behind the step size adaptation algorithm? -- Also, why is it that you increase x before the step calculations? – LutzL May 19 '17 at 9:25
• @LutzL the equation for step size was taken from a professor's lecture I took in my 3rd year undergrad numerical methods course. As for increasing the x value, I believe that was a programming logic error. – Curtis May 19 '17 at 22:24
• There should be some more background, like y5-y4 being $Ch^5$ with $C$ independent of $h$ and slowly changing with $x$, that the local estimate for the global error is then $Ch^4$ which should be close to the value delta0 from the relative error tolerance etc. This should then show how the used formulas and logic have to be corrected. Also, there is no provision in the logic to increase the step size if the error estimate is too small. – LutzL May 20 '17 at 7:27

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.

• Thank you very much for your very thorough answer! Would you mind expanding a bit on your explanation for the optimal step size formula? – Curtis May 19 '17 at 22:07