|
#! /usr/bin/env python
# -*- coding: utf-8 -*-
"""
Examples of several algorithm for numerical integration problems.
See the `integrals.py <integrals.html>`_ file for more details.
- *Date:* Saturday 18 June 2016, 18:59:23.
- *Author:* `Lilian Besson <https://bitbucket.org/lbesson/>`_ for the `CS101 course <http://perso.crans.org/besson/cs101/>`_ at `Mahindra Ecole Centrale <http://www.mahindraecolecentrale.edu.in/>`_, 2015.
- *Licence:* `MIT Licence <http://lbesson.mit-license.org>`_, © Lilian Besson.
"""
from __future__ import division, print_function, absolute_import # Python 2/3 compatibility
import math
from integrals import *
# Do some quick tests (it should be in the 'tests.py' file):
19 ↛ exitline 19 didn't exit the module, because the condition on line 19 was never falseif __name__ == '__main__':
print("\nSome examples of numerical integrations:\n")
# # %% Example 1
f = math.cos
print("Plotting and computing the integral of cos(x) from 0 to pi.")
print("Formally, the value is about 0 (= sin(pi) - sin(0)).")
print("With 15 left rectangles:", riemann_left(f, 0, math.pi, n=15))
print("With 15 right rectangles:", riemann_right(f, 0, math.pi, n=15))
print("With 15 center rectangles:", riemann_center(f, 0, math.pi, n=15))
print("With 3 trapezoides:", trapez(f, 0, math.pi, n=3))
# print("With 15 left rectangles:", plot_riemann_left(f, 0, math.pi,
# namef=r"$x \mapsto \cos(x)$",n=15, figname="riemannleft.png"))
# print("With 15 right rectangles:", plot_riemann_right(f, 0, math.pi,
# namef=r"$x \mapsto \cos(x)$",n=15, figname="riemannright.png"))
# print("With 15 center rectangles:", plot_riemann_center(f, 0, math.pi,
# namef=r"$x \mapsto \cos(x)$",n=15, figname="riemanncenter.png"))
# print("With 3 trapezoides:", plot_trapez(f, 0, math.pi,
# namef=r"$x \mapsto \cos(x)$",n=3, figname="trapezoides.png"))
# # %% Example 2
# f = lambda x: x**2
# print("Plotting and computing the integral of x**2 from 0 to 3.")
# print("Formally, the value is about 9 (= 3**3 / 3).")
# print("With 8 left rectangles:", plot_riemann_left(f, 0, 3,
# namef="$x \mapsto x^2$", n=8, figname="riemannleft2.png"))
# print("With 8 right rectangles:", plot_riemann_right(f, 0, 3,
# namef="$x \mapsto x^2$", n=8, figname="riemannright2.png"))
# print("With 8 center rectangles:", plot_riemann_center(f, 0, 3,
# namef="$x \mapsto x^2$", n=8, figname="riemanncenter2.png"))
# print("With 3 trapezoides:", plot_trapez(f, 0, 3,
# namef="$x \mapsto x^2$", n=3, figname="trapezoides2.png"))
# %% Example 3 with Monte-Carlo
n = 10000
xmin, xmax = 1, 6
def f1(x):
""" f1(x) = x"""
return x
# For this particular example, its easy:
ymin, ymax = xmin, xmax
print("\nFor this function f1: x → x, on I = [{}, {}], we will pick {} random points.".format(xmin, xmax, n))
print("Manually, I chose ymin = {}, ymax = {}.".format(ymin, ymax))
random.seed(1)
intf = montecarlo(f1, xmin, xmax, n, ymin, ymax)
print("This leads to an approximated integral value of {}.".format(intf))
print("Error % is", 100 * abs(17.5 - intf))
ymin, ymax = yminmax(f1, xmin, xmax, n) # guess!
print("But experimentally, I found a possible ymin, ymax to be {}, {}.".format(ymin, ymax))
intf = montecarlo(f1, xmin, xmax, n, ymin, ymax)
print("This leads to a second approximated integral value of {}.".format(intf))
print("Error % is", 100 * abs(17.5 - intf))
# plot_montecarlo(f, xmin, xmax, 1500, ymin, ymax, namef=r"$x \mapsto x$", figname="montecarlo.png")
print("Formally, we compute the integral as 17.5 (36/2 - 1/2).")
# %% Example 4 with Monte-Carlo
xmin, xmax = 0.0, 1.0
def f2(x):
""" f1(x) = x"""
return x ** 3
# For this particular example, its easy:
ymin, ymax = xmin, xmax
print("\nFor this function f2: x → x^3, on I = [{}, {}], we will pick {} random points.".format(xmin, xmax, n))
print("Manually, I chose ymin = {}, ymax = {}.".format(ymin, ymax))
intf = montecarlo(f2, xmin, xmax, n, ymin, ymax)
print("This leads to an approximated integral value of {}.".format(intf))
print("Error % is", 100 * abs(0.25 - intf))
ymin, ymax = yminmax(f2, xmin, xmax, n) # guess!
print("But experimentally, I found a possible ymin, ymax to be {}, {}.".format(ymin, ymax))
intf = montecarlo(f2, xmin, xmax, n, ymin, ymax)
print("This leads to a second approximated integral value of {}.".format(intf))
print("Error % is", 100 * abs(0.25 - intf))
# plot_montecarlo(f2, xmin, xmax, ymin, ymax, 1500, namef=r"$x \mapsto x^3$", figname="montecarlo2.png")
print("Formally, we compute the integral as 0.25 (1/4).")
# %% Example 5 with Monte-Carlo
xmin, xmax = 1e-8, 3.
def f3(x):
""" f3(x) = x"""
return 1.0 / (1 + (math.sinh(2 * x)) * (math.log(x)**2))
# For this particular example, its easy:
ymin, ymax = 0, 1
print("\nFor this function f3: x → 1/(1+sinh(2x)log(x)^2) on I = [{}, {}], we will pick {} random points.".format(xmin, xmax, n))
print("Manually, I chose ymin = {}, ymax = {}.".format(ymin, ymax))
intf = montecarlo(f3, xmin, xmax, n, ymin, ymax)
print("This leads to an approximated integral value of {}.".format(intf))
ymin, ymax = yminmax(f3, xmin, xmax, n) # guess!
print("But experimentally, I found a possible ymin, ymax to be {}, {}.".format(ymin, ymax))
intf = montecarlo(f3, xmin, xmax, n, ymin, ymax)
print("This leads to a second approximated integral value of {}.".format(intf))
# plot_montecarlo(f3, xmin, xmax, ymin, ymax, 1500, namef=r"$x \mapsto 1/(1+\sinh(2x)\log(x)^2)$", figname="montecarlo3.png")
# plot_montecarlo(f3, xmin, xmax, ymin, ymax, 10000, namef=r"$x \mapsto 1/(1+\sinh(2x)\log(x)^2)$", figname="montecarlo4.png")
# End of tests.py
|
