Intégration numérique

Voici quelques éléments pour vous aidez à répondre aux question du TP6.

Les autres TP sont à retrouver à la page d’accueil

Ce TP est l’occasion d’utiliser les modules vu à l’exercice 5 du TP4. Les modules servent à structurer votre code. Rappel: Je ne cherche pas à corriger le TP6 ici. Ainsi toutes les questions ne seront pas corrigées mais vous trouverez un programme vous permettant de répondre par vous même à toutes les questions. Réfléchissez, piochez avec sagesse, adaptez !

Méthode de quadrature

module calc_int
    use fonction

    implicit none
contains

    subroutine set_method(method, filename, W, T)
        integer               , intent(in   ) :: method
        character(16)         , intent(  out) :: filename
        real(PR), dimension(:), allocatable, intent (inout) :: W, T

        select case(method)

        case (0) !"trapeze"
            filename = "erreur_trapz.dat"
            w = (/1.,1./)
            T = (/-1.,1./)

        case (1) ! "middle"
            filename = "erreur_ptmil.dat"
            W = (/2._pr/)
            T = (/0._pr/)

        case (2) ! "Simpson"
            filename = "erreur_simps.dat"
            W = (/ 1._pr, 4._pr, 1._pr/)
            W = W/3._pr
            T = (/-1._pr, 0._pr, 1._pr/)

        ! Gauss-Legendre
        case (3) ! Gauss 2 pts
            filename = "erreur_ga2pt.dat"
            W = (/1._pr, 1._pr/)
            T = (/-sqrt(3._pr)/3._pr, sqrt(3._pr)/3._pr/)

        case (4) ! Gauss 3 pts
            filename = "erreur_ga3pt.dat"
            W = (/5._pr, 8._pr, 5._pr/)
            W = W/9._pr
            T = (/-sqrt(3._pr/5), 0._pr, sqrt(3._pr/5) /)

        ! Gauss-Laguerre
        case (5) ! Gauss-Laguere 1 pts
            W = (/2._pr/)
            T = (/1._pr/)

        case (6) ! Gauss-Laguere 2 pts
            W = (/2._pr+sqrt(2._pr), 2-sqrt(2._pr)/)
            W = W/4._pr
            T = (/2._pr-sqrt(2._pr), 2+sqrt(2._pr)/)

        case default ! rectangle a gauche
            filename = "erreur_rectg.dat"
            W = (/2._pr/)
            T = (/-1._pr/)
        end select

    end subroutine set_method

    function quadrature(a, b, n, W, T) result(approx)
        real(PR)              , intent(in) :: a, b
        integer(PR)           , intent(in) :: n
        real(PR), dimension(:), intent(in) :: W, T
        ! returns
        real(PR) :: approx

        ! locals variables
        real(PR) :: ai, h
        integer  :: i, l


        h = (b-a)/n
        approx = 0.
        do i = 0, n-1
            ai = a + i*h
            do l=1, size(T)
                approx = approx + W(l) * f( ai + h/2.*(1+T(l)) )
            end do
        end do
        approx = approx * h/2
    end function quadrature

end module calc_int
Dans ce premier module, je rassemble toutes les fonctions qui permettent de définir la quadrature.

set_method renvoit W et T (en les allouants) mais aussi filename, qui me servira à faire mes graphiques avec gnuplot.

quadrature renvoit la quadrature

Programme principal

program part2
    use calc_int
implicit none

    real(PR), parameter :: pi = 4*atan(1.)

    real(PR), dimension(:), allocatable :: W, T

    real(PR)            :: approx
    integer             :: method, n
    real(PR) :: a, b
    ! Donnee exacte
    real(PR) :: ExactValue

    character(16) :: filename

    print *, "Entrez la borne min :"
    read *, a
    print *, "Entrez la borne max :"
    read *, b
    ! print *, "Entrez le nombre de points :"
    ! read *, nmax

    print *, "Choisir une méthode:"
    print *, "   0 : trapeze"
    print *, "   1 : milieu"
    print *, "   2 : Simpson"
    print *, "   3 : Gauss 2pts"
    print *, "   4 : Gauss 3pts"
    print *, "   autre: rectangle gauche"
    read *, method

    ExactValue = integral(a, b)

    ! W et T seront alloues par set_method
    call set_method(method, filename, W, T)
    open(11, file=filename)

    do n=0, 15
        approx = quadrature(a, b, 2_pr**n, W, T)
        print "(i3, 1x, i15, E15.7, 1x, i2)", n, 2_pr**n, abs(approx-ExactValue), method
        ! dans le fichier on retrouve
        ! colonne 1: h, colonne 2: l'erreur, colonne 3: le nombre d'évaluation de f
        write(11, "(2E15.7, 1x, i15)")  (b-a)/(2**n), abs(approx - ExactValue), 2**n * size(T)
    end do

    close(11)
    deallocate(W, T)

end program part2
! set logscale xy
! set key bottom
! plot 'erreur_rectg.dat' w lp, 'erreur_ptmil.dat' w lp, 'erreur_trapz.dat' w lp, 'erreur_simps.dat' w lp, 'erreur_ga2pt.dat' w lp, 'erreur_ga3pt.dat' w lp
!plot 'erreur_rectg.dat' using 1:2 w lp, 'erreur_ptmil.dat' using 1:2 w lp, 'erreur_trapz.dat' using 1:2 w lp, 'erreur_simps.dat' using 1:2 w lp, 'erreur_ga2pt.dat' using 1:2 w lp, 'erreur_ga3pt.dat' using 1:2 w lp
!plot 'erreur_rectg.dat' using 3:2 w lp, 'erreur_ptmil.dat' using 3:2 w lp, 'erreur_trapz.dat' using 3:2 w lp, 'erreur_simps.dat' using 3:2 w lp, 'erreur_ga2pt.dat' using 3:2 w lp, 'erreur_ga3pt.dat' using 3:2 w lp

fonction

module fonction
implicit none
    integer, parameter :: PR = 8

contains
    function f(t) result(rslt)
        real(PR), intent(in) :: t
        ! Locals
        real(PR) :: rslt

        rslt = 1._pr/t
    end function f


    function exact(t) result(rslt)
        real(PR), intent(in) :: t
        ! returns
        real(PR) :: rslt

        ! rslt = log(t)
        rslt = log(t)
    end function exact

    function integral(a, b) result(rslt)
        real(PR), intent(in) :: a, b
        ! returns
        real(PR) :: rslt

        rslt = exact(b)-exact(a)
    end function integral
end module fonction
Le module le plus basique. Il contient tout ce qui est lié à la fonction f. Ici f(x)=1/x. On retrouve bien une organisation par thématique, simple, épurée et minimaliste. Cela apporte de la lisibilité.

Comment compiler

$ gfortran fonction.f90 calc_int.f90 main2.f90 -o main2.exe
$ ./main2.exe

Tracer l’évolution de l’erreur avec Gnuplot

Rappelons que mes fichiers .dat contiennent 3 colonnes. Il font donc préciser à Gnuplot celles qu’il devra utiliser pour les tracer

set title "Evolution erreur absolue"
set xlabel "pas"
set ylabel "erreur"
set key bottom
plot 'erreur_rectg.dat' using 1:2 w lp, 'erreur_ptmil.dat' using 1:2 w lp, 'erreur_trapz.dat' using 1:2 w lp, 'erreur_simps.dat' using 1:2 w lp, 'erreur_ga2pt.dat' using 1:2 w lp, 'erreur_ga3pt.dat' using 1:2 w lp

La conclusion est évidente et tient sur 2 points (Meilleure approximation avec un pas de temps diminuant jusqu’à un certain point. Pour certaines méthodes, nous voyons que passer un certain cap l’erreur reste faible mais augmente)

set title "Evolution erreur absolue"
set xlabel "nb evaluation f"
set ylabel "erreur"
set key left bottom
plot 'erreur_rectg.dat' using 3:2 w lp, 'erreur_ptmil.dat' using 3:2 w lp, 'erreur_trapz.dat' using 3:2 w lp, 'erreur_simps.dat' using 3:2 w lp, 'erreur_ga2pt.dat' using 3:2 w lp, 'erreur_ga3pt.dat' using 3:2 w lp

Déterminer le nombre d’itérations pour atteindre une erreur spécifique

En utilisant ce que nous avons déjà, il est possible d’itérer jusqu’à obtenir l’erreur acceptable. Mais ne partons pas à n = 0

Nous pouvons modifier une partie du code (allez-vous la trouver ?) par celle-ci (à compléter):

n = nbeg
do while ((erreur > stopCriteria) .and. (n<29))
    ...
    n = n+1
end do

Deux nouvelles variables sont à demander à l’utilisateur

  • nbeg est la puissance de 2 qui détermine où commencer la boucle,
  • stopCriteria est l’erreur acceptable.

Pour déterminer nbeg, il suffit de regarder la courbe de l’erreur en fonction du pas h et de la suivre voir vidéo en dessous.

On obtient le tableau suivant (attention le nombre de subdivision est 2**n et non n)

1,00E-06 1,00E-09 1,00E-12
h< n h< n h< n
rectangle0,001 6 ? ? ? ?
trapz 0,0042065 0,001 6 ? ?
ptmil 0,00667 5 0,001 6 ? ?
simpson 0,1385 1 0,03134 3 0,004905 5
gauss 2pt0,1385 1 0,03134 3 0,00490 5
gauss 3pt0,5524 0 0,1826 1 0,05451 2

Certaines cases ne peuvent être remplies uniquement avec le graphique.