Как перебрать цикл по epsilon = 1e-8 для реализации интегратора Симпсона - PullRequest
Новая
       68

Как перебрать цикл по epsilon = 1e-8 для реализации интегратора Симпсона

0 голосов
/ 19 октября 2019

Я реализовал следующую логику и задал этот вопрос для другого вопроса (диапазон массива). Я получаю вывод, но он не проходит цикл for для итерации, потому что я дал frange(start, stop, range)

Объяснение 

"""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
"""

Код : 

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
%matplotlib inline
import pylab as pl
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  = 2 # 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)

    while error > epsilon:
        for i in pl.frange(1,epsilon,1):
            print('Hello')
            if(i%2 ==0):
                print('Hello')
                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)

Другая функция, вызывающая вышеупомянутую функцию 

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) )

Ниже приводится результат 

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

Это должно дать мне больше итераций

Может кто-нибудь помочь мне перебрать цикл, чтобы получить больше результата

...