[2022 Fall] Final Exam - Example Code¶

This is the example code from TAs.

Romberg Intergration Example¶

In [1]:

from scipy import integrate
from scipy.special import erf
import numpy as np
gaussian = lambda x: ((2*x) + (3/x))**2
result = integrate.romberg(gaussian, 1, 2, show=True)

from scipy import integrate from scipy.special import erf import numpy as np gaussian = lambda x: ((2*x) + (3/x))**2 result = integrate.romberg(gaussian, 1, 2, show=True)

Romberg integration of <function vectorize1.<locals>.vfunc at 0x7fe76cbdbb50> from [1, 2]

 Steps  StepSize   Results
     1  1.000000 27.625000 
     2  0.500000 26.312500 25.875000 
     4  0.250000 25.955944 25.837092 25.834565 
     8  0.125000 25.864188 25.833602 25.833370 25.833351 
    16  0.062500 25.841060 25.833351 25.833334 25.833333 25.833333 
    32  0.031250 25.835266 25.833334 25.833333 25.833333 25.833333 25.833333 

The final result is 25.833333333536203 after 33 function evaluations.

Trapeziodal Rule Example¶

In [2]:

x=np.array([0,1,2,3])
fx=(x**2)*(np.exp(x))
np.trapz(fx, dx=1)

x=np.array([0,1,2,3]) fx=(x**2)*(np.exp(x)) np.trapz(fx, dx=1)

Out[2]:

122.65942237852616

Simpson Rule Example (composite simpson (1/3))¶

In [3]:

integrate.simps(fx,dx=1)

integrate.simps(fx,dx=1)

Out[3]:

108.27484014325017

Gauss Quadrature Example¶

In [4]:

from scipy import integrate

f = lambda x: 1/(1+x**2)
integrate.quadrature(f, -3.0, 3.0)

from scipy import integrate f = lambda x: 1/(1+x**2) integrate.quadrature(f, -3.0, 3.0)

Out[4]:

(2.498091551681014, 2.0135778822094608e-08)

Richardson Extropolation¶

In [5]:

import numpy as np

def f(x): 
    return np.cos(x)

def fp(x): 
    return np.sin(x)

fp(np.pi/4)

def phi(x,h): 
    return (f(x+h)-f(x-h))/(2*h)

# d = [phi(1,h) for h in [2**(-n) for n in range(5)]]
d = [phi(np.pi/4,h) for h in [np.pi/3,np.pi/6]]

import numpy as np def f(x): return np.cos(x) def fp(x): return np.sin(x) fp(np.pi/4) def phi(x,h): return (f(x+h)-f(x-h))/(2*h) # d = [phi(1,h) for h in [2**(-n) for n in range(5)]] d = [phi(np.pi/4,h) for h in [np.pi/3,np.pi/6]]

In [6]:

d

d

Out[6]:

[-0.5847726009252571, -0.6752372371178295]

In [7]:

N = len(d)
D = np.zeros((N,N))
D[:,0] = d
for m in range(1,N):
    for n in range(m,N):
        D[n,m] = (4**m*D[n,m-1]-D[n-1,m-1])/(4**m-1)

N = len(d) D = np.zeros((N,N)) D[:,0] = d for m in range(1,N): for n in range(m,N): D[n,m] = (4**m*D[n,m-1]-D[n-1,m-1])/(4**m-1)

In [8]:

D

D

Out[8]:

array([[-0.5847726 ,  0.        ],
       [-0.67523724, -0.70539212]])

---

Last update: 2024-04-27