lab10-extra

trapez(f, a, b, n)

Composite trapezoidal integration

Background

The equation below shows that the integral \(I\) of the function \(f(x)\) from \(x=a\) to \(x=b\) can be approximated by \(A\) through the so-called composite trapezoidal integration rule:

\[I = \int\limits_a^b f(x) \mathrm{d}x \approx A = \frac{h}{2}\left(f(a) + f(b) + 2\sum\limits_{i=1}^{n-1}f(x_i)\right)\]

where \(h=\frac{b-a}{n}\), and \(x_i=a+ih\), and \(n\) is the number of subdivisions.

Example:

Composite trapezoidal rule with 4 subdivisions.

The figure demonstrates the idea for \(f(x)=x^2\), \(a=0\), \(b=2\), \(n=4\) and thus \(h=0.5\):

The composite trapezoidal rule computes the area under the red line exactly, where the red line is an approximation of the actual function of interest (here \(f(x)=x^2\)) shown by a black dashed line.

Exercise (trapez):

Write a function trapez(f, a, b, n) which is given

  • a function f which depends on one parameter (i.e. f(x)) and returns a number

  • a lower (a) and upper (b) integration limit

  • and the number of subdivisions (n) to be used.

The function should use the composite trapezoidal rule to compute A and return this value.

Examples

In [ ]: def f(x):
   ...:     return x
   ...:

In [ ]: trapez(f, 0, 1, 1)
Out[ ]: 0.5

In [ ]: trapez(f, 0, 1, 5)
Out[ ]: 0.5

In [ ]: def f2(x):
   ...:     return x * x
   ...:

In [ ]: trapez(f2, 0, 1, 1)
Out[ ]: 0.5

In [ ]: trapez(f2, 0, 1, 10)
Out[ ]: 0.3350000000000001

In [ ]: trapez(f2, 0, 1, 100)
Out[ ]: 0.33335000000000004

In [ ]: trapez(f2, 0, 1, 1000)
Out[ ]: 0.33333349999999995

finderror(n)

Estimate the numerical integration error

Write a function finderror(n) which uses the trapez() function to find the error of the trapezoidal integral approximation. The function should compute the integral of \(f(x) = x^2\) with integration limits \(a = 0\) and \(b = 2\) numerically. The function should then subtract this numerical result \(A\) from the exact integral value \(I\) and return the difference. Using anytical methods (i.e. symbolic integration), we can compute the exact integral to be \(I=\int\limits_0^2f(x)\mathrm{d}x = \frac{8}{3}\).

Example:

In [ ]: finderror(5)
Out[ ]: -0.05333333333333412

(You should find approximately that the error decays by a factor 4 if you double n (because the trapez method has an error of order \(h^2\) and \(h\) is proportional to \(1/n\).)

using_quad()

Numerical integration

Write a function using_quad() that computes and returns the integral \(\int\limits_a^b x^2 f(x) \mathrm{d}x\) with \(f(x) = x^2\) from \(a=0\) to \(b=2\) using scipy.integrate.quad. The function using_quad should return the value that the quad() function returns, i.e. a tuple (y, abserr) where y is quad’s best approximation of the true integral, and abserr an absolute error estimate.

You should find that this method (with the default settings) is (far) more accurate than our trapez function.

Please include the extra tasks in your file lab10.py and submit as lab10 assignment.

Back to lab10.

End of lab10-extra.