# -*- coding: utf-8 -*-
"""
Created on Wed Jan 07 16:20:59 2015

@author: grivet

Listing 8.1 : intégration de Romberg
Le module 'scipy.integrate' comporte un programme 'romberg'
"""
import numpy as np

def affiche_mat(titre,M):
    print (titre)
    for row in M:
        for val in row:
            print( ' %10.5f' %val, end = '')
        print()

def fn(x):
    return np.exp(x)*np.cos(x)

lmax = 5
a = 0; b = np.pi
h = b-a
J = np.zeros((lmax,lmax))
J[0,0] = 0.5*h*(fn(a)+fn(b))
for l in range(1,lmax):
    h = h/2
    x = np.arange(a+h,b,2*h)
    di = h*sum(fn(x))
#    print di
    J[l,0] = 0.5*J[l-1,0] + di
for c in range(1,lmax):
    for l in range(c,lmax):
        J[l,c] = (4**(c)*J[l,c-1]-J[l-1,c-1])/(4**(c)-1)
print ('J exacte = %12.7f' % (-(np.e**np.pi+1)/2) )
affiche_mat('Romberg: ',J)
print ('J calculée = %12.7f' % J[-1,-1])