0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – LutzL May 19 '17 at 9:25
  • $\begingroup$ @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. $\endgroup$ – Curtis May 19 '17 at 22:24
  • $\begingroup$ 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. $\endgroup$ – LutzL May 20 '17 at 7:27
2
$\begingroup$

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.

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

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.