5.5. Exercise solutions

3.2.1

import numpy as np
import scipy.integrate as integrate
N = 1000000
L = 10000.
h = L/N
x = np.linspace(1,L,N+1)
f = 1./x**2
integrate.trapezoid(f,x=x)

0.9999166630003236

3.2.2

Using the transformation \(x=(1-y)/y\),

\[ \int_0^\infty \frac{1}{x^2+1} \, dx = \int_0^1 \frac{1}{2y^2-2y+1}\, dy \]

Evalute the right hand side using the Simpson’s rule with \(h=0.01\).

import numpy as np
import scipy.integrate as integrate

N=100
h=1./N

y = np.linspace(0,1,N+1)
f = 1.0/(2*y**2 - 2*y + 1)
integral=integrate.simpson(f,x=y)
error=integral-np.pi/2
print("simpson =",integral)
print("  exact =",np.pi/2.)
print("  error =",error)

simpson = 1.570796326793627
  exact = 1.5707963267948966
  error = -1.269651050961329e-12

Notice that the error is slightly bigger than the result of Example 3.2.1.

3.2.3

import numpy as np
import scipy.integrate as integrate

def f(x):
    return np.sin(x)*np.exp(-x)

gauss=integrate.quad(f,0,np.inf)

print(gauss)

(0.5000000000000002, 1.4875911973447031e-08)

import numpy as np

# We use 4-point method.
x=np.array([3.2254768961939231180e-1, 1.7457611011583465757,
            4.5366202969211279833,    9.3950709123011331292])

w=np.array([6.0315410434163360164e-1, 3.5741869243779968664e-1,
            3.8887908515005384272e-2, 5.3929470556132745010e-4])    

gauss=(w*np.sin(x)).sum()
exact=1/2.

print(" Exact={0:18.12e}\n Gauss={1:18.12e}\n Error={2:18.12e}".format(exact,  gauss, abs(exact-gauss)))

 Exact=5.000000000000e-01
 Gauss=5.048792794602e-01
 Error=4.879279460199e-03

We got only 2 digits of precision. It is still remarkable considering the function is evaluated only at four points over infinite interval.