0

I have tried to implement this question in my logic to calculate Simpson's Method Adaptive Implementation using epsilon = 1e-8. Following is the explanation:

"""The approximate definite integral of function from a to b using Simpson's method. 
This function is vectorized, it uses numpy array operations to calculate the approximation.
This is an adaptive implementation, the method starts out with N=2 intervals, and try
successive sizes of N (by doubling the size), until the desired precision, is reached.
This adaptive solution uses our improved approach/equation for Simpson's method, to avoid 
unnecessary recalculations of the integrand function.

a, b - Scalar float values, the begin, and endpoints of the interval we are to
        integrate the function over.
f - A vectorized function, should accept a numpy array of x values, and compute the
      corresponding y values for all points according to some function.
epsilon - The desired precision to calculate the integral to.  Default is 8 decimal places
      of precision (1e-8)

returns - A tuple, (ival, error).  A scalar float value, the approximated integral of 
       the function over the given interval, and a scaler float value of the 
       approximation error on the integral"""

Following is my code for Simpson's Method:

import pylab as pl
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
%matplotlib inline
def simpsons_adaptive_approximation(a, b, f, epsilon=1e-8):

    N_prev = 2 # the previous number of slices
    h_prev = (b - a) / N_prev # previous interval width
    x = np.arange(a+h_prev, b, h_prev) # x locations of the previous interval
    I_prev =  h_prev * (0.5 * f(a) + 0.5 * f(b) + np.sum(f(x)))

    # set up variables to adaptively iterate successively better approximations
    N_cur  = 4 # the current number of slices
    I_cur  = 0.0 # calculated in loop iteration
    error  = 1.0 # calculated in loop iteration
    itr    = 1 # keep track of the number of iterations we perform, for display/debug

    h = (b-a)/float(epsilon)
    I_cur = f(a) + f(b)

    for i in frange(1,epsilon,1):
        if(i%2 ==0):
            I_cur = I_cur + (2*(f(a + i*h)))
        else:
            I_cur = I_cur + (4*(f(a + i*h)))
    error = np.abs((1.0/3.0) * (I_cur - I_prev))
    print("At iteration %d (N=%d), val=%0.16f prev=%0.16f error=%e" % (itr, N_cur, I_cur, I_prev, error) )

    I_cur *= (h/3.0)

    I_prev = I_cur
    N_prev = N_cur
    N_cur *= 2
    itr += 1

return (I_cur, error)

Following is my f2(x) function:

def f2(x):
    return x**4 - 2*x + 1

a = 0.0
b = 2.0
eps = 1e-10
(val, err) = simpsons_adaptive_approximation(a, b, f2, eps)
print( "Calculated value: %0.16f  error: %e  for an epsilon of: %e" % (val, err, eps) )

**outcome is only 1 iteration **:

At iteration 1 (N=4), val=14.0000000000000000 prev=7.0000000000000000 error=2.333333e+00
Calculated value: 93333333333.3333435058593750  error: 2.333333e+00  for an epsilon of: 1.000000e-10

And Expected Result Something like this for multiple iteration:

At iteration 1 (N=...), val=........ prev=...... error=1.552204e-10
At iteration n (N=...), val=........ prev=...... error=3.880511e-11
Calculated value: 0.0000000000  error: 3.880511e-11  for an epsilon of: 1.000000e-10

I want to explore the code with some understanding because this question is already asked for other values as n but not for an array

Mayur Potdar
  • 451
  • 1
  • 8
  • 18
  • Possible duplicate of [range() for floats](https://stackoverflow.com/questions/7267226/range-for-floats) – walnut Oct 19 '19 at 00:24
  • @uneven_mark I took a look at your suggestion but It's not satisfying my understanding for the above code – Mayur Potdar Oct 19 '19 at 00:27
  • `range` does only accept integer values. I don't know what more there is to explain the error. The duplicate and links in the duplicate provide alternative approaches for `range` with floats. The rest of the code you posted does not relate to the error as far as I can tell. – walnut Oct 19 '19 at 00:31
  • @uneven_mark Hey, I was looking at multiple solution to understand and I fixed my code. Now I'm getting following result `At iteration 1 (N=4), val=14.0000000000000000 prev=7.0000000000000000 error=2.333333e+00 Calculated value: 93333333333.3333435058593750 error: 2.333333e+00 for an epsilon of: 1.000000e-10` I should have obtain more results before calculated value. – Mayur Potdar Oct 19 '19 at 00:33
  • @uneven_mark Somehow It's not going in for loop and just doing final iteration can you please help me to understand? – Mayur Potdar Oct 19 '19 at 00:54
  • Please don't remove/replace the erroneous code in your question. It makes it impossible for other readers to follow. Also, if you have a new question unrelated to the original one, please post a new question. In this particular case though, note that `frange(1,1e-10,1)` does not do any iterations, because `start (1) < end (1e-10)`. – walnut Oct 19 '19 at 01:00
  • @uneven_mark I apologize and I was trying to neglect the confusion in the code. Can you please tell me how can I iterate over each value to calculate Simpson's adaptive implementation – Mayur Potdar Oct 19 '19 at 01:08

0 Answers0