rombergSoluce.py

# -------------------------------------------------------------------------
#
# PYTHON for DUMMIES 19-20
# Problème 4
#
# Solution détaillée de notre implémentation de Romberg :-)
#  Vincent Legat
#
# -------------------------------------------------------------------------
# 

from numpy import *

def romberg(f,a,b,n,nmax,tol) :

  size = int(ceil(log2(nmax/n)))
  I = zeros((size,size))
  
  i = j = 0
  X,h = linspace(a,b,n+1,retstep=True)
  I[i,j] = trapz(f(X),X) 
  print('\n      Step %2d %2d : I = %21.14e  (n = %d) ' % (i,j,I[i,j],n))

  errorEst = tol + 1; i = 0
  while n*2 < nmax and errorEst > tol: 
    i = i+1; j = 0; hold = h; n = n*2; h = h/2
    X = arange(a+h,b,hold) 
    I[i,j] = (I[i-1,j]/hold + sum(f(X)))*h
    print('\n      Step %2d %2d : I = %21.14e  (n = %d)' % (i,j,I[i,j],n))
    for j in range(1,i+1):
      coeff = 2**(2*j)
      I[i, j] = (coeff * I[i,j-1] - I[i-1,j-1])/(coeff - 1)      
      print('      Step %2d %2d : I = %21.14e' % (i,j,I[i,j]))
    errorEst = abs(I[i,j] - I[i,j-1])  
    
  return I[i,j],n,errorEst