#! /usr/bin/env python # -*- coding: utf-8 -*- """ Examples of several algorithm for numerical integration problems. See the `integrals.py `_ file for more details. :Date: Tuesday, April 28, 08:45:29 2015. :Author: `Lilian Besson `_ for the `CS101 course `_ at `Mahindra Ecole Centrale `_, 2015. :Licence: `MIT Licence `_, © Lilian Besson. """ from __future__ import division import math from integrals import * # Do some quick tests (it should be in the 'tests.py' file): if __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 f(x): return x # f = lambda x: x # quicker! # For this particular example, its easy: ymin, ymax = xmin, xmax print "\nFor this function f: 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(f, 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(f, xmin, xmax, n) # guess! print "But experimentally, I found a possible ymin, ymax to be {}, {}.".format(ymin, ymax) intf = montecarlo(f, 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., 1. f = lambda x: x**3 # quicker! # For this particular example, its easy: ymin, ymax = xmin, xmax print "\nFor this function f: 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(f, 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(f, xmin, xmax, n) # guess! print "But experimentally, I found a possible ymin, ymax to be {}, {}.".format(ymin, ymax) intf = montecarlo(f, 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(f, 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. f = lambda x: 1./(1+(math.sinh(2*x))*(math.log(x)**2)) # quicker! # For this particular example, its easy: ymin, ymax = 0, 1 print "\nFor this function f: 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(f, xmin, xmax, n, ymin, ymax) print "This leads to an approximated integral value of {}.".format(intf) ymin, ymax = yminmax(f, xmin, xmax, n) # guess! print "But experimentally, I found a possible ymin, ymax to be {}, {}.".format(ymin, ymax) intf = montecarlo(f, xmin, xmax, n, ymin, ymax) print "This leads to a second approximated integral value of {}.".format(intf) # plot_montecarlo(f, xmin, xmax, ymin, ymax, 1500, namef=r"$x \mapsto 1/(1+\sinh(2x)\log(x)^2)$", figname="montecarlo3.png") # plot_montecarlo(f, xmin, xmax, ymin, ymax, 10000, namef=r"$x \mapsto 1/(1+\sinh(2x)\log(x)^2)$", figname="montecarlo4.png") # End of tests.py