Simpson revival

[Namir]

(10-02-2016 03:29 PM)Dieter Wrote:  

(09-30-2016 06:00 AM)I Wrote:  The test should compare the last and previous improved approximations, cf. the posted VBA code. In your program that's the current and previous value of Y.

It looks like even this remaining error can be estimated and an even more accurate result can be calculated from the last two improved values. Here is a more structured experimental version.

Code:

Function Simpson(a, b, Optional errlimit = 0.000001, Optional nmax = 3000)
   
   fab = f(a) + f(b)
   s2 = 0
   s4 = 0
   simp_new = 0
   improved_new = 0
   n = 2
Do
   h = (b - a) / n
   s2 = s2 + s4
   s4 = 0
   
   For i = 1 To n - 1 Step 2
      s4 = s4 + f(a + i * h)
   Next
   
   simp_old = simp_new
   improved_old = improved_new
   simp_new = (fab + 2 * s2 + 4 * s4) * h / 3
   
   If n = 2 Then
      improved_new = simp_new
      converged = False
   Else
      improved_new = (16 * simp_new - simp_old) / 15
      converged = Abs(improved_new - improved_old) < 63 * errlimit
   End If
   
   n = 2 * n

Loop Until converged Or n > nmax

Simpson = (64 * improved_new - improved_old) / 63

End Function

Function f(x)
   f = 1 / x
End Function

What do you think?

The following diagram may give an impression for the integral of f(x) = 1/x from a=1 to b=2. Here X is the number of iterations (i.e. n=2^x intervals) and Y is the accuracy, i.e. the number of valid digits. The red line is the classic Simpson method, the blue one represents the "improved" method, and finally the green graph shows the further extrapolated result as shown in the last program.



Dieter

You latest code works in general better than Pekis' version. The ratio of errors generated by your code to that of Pekis range in the orders [1, 10]. Good work!