Le programme

Voici le programme FORTRAN correspondant à la méthode de Newton pour résoudre le système non linéaire (8.1). Nous avons choisi $f$ la fonction définie par $c(x)=sin(x)$. La solution dans ce cas précis est donnée par $x_1=x_2=\cdots=x_N=0$.

        program newton
        implicit double precision (a-h,o-z)
        parameter (nmax=1000)
        dimension x(nmax),dl(nmax-1),d(nmax),du(nmax-1),b(nmax)
c
c       initialisations
c
        n=10
        do 10 i=1,n
                x(i) = 1.
10        continue
c
c       methode de Newton : 
c          resolution du syst. lin. : Df( x^n ) ( y^n ) = -f( x^n )
c          correction : x^n+1 = x^n + y^n
c
        itmax=20
        print*,' iterations'    
        do 40 iter=1,itmax
c
c          construction de la matrice
c
           dl(1) =-1.
           d(1)  =2.+c(x(1),1)
           du(1) =-1.
           b(1)  =-2.*x(1)+x(2)-c(x(1),0)
           do 20 i=2,n-1
                dl(i) =-1.
                d(i)  =2.+c(x(i),1)
                du(i) =-1.
                b(i)  =x(i-1)-2.*x(i)+x(i+1)-c(x(i),0)
20         continue
           d(n) =2.+c(x(n),1)
           b(n) =-2.*x(n)+x(n-1)-c(x(n),0)
c
c          resolution du syst. lin. : appel de DGTSV (LAPACK)
c 
           call dgtsv(n,1,dl,d,du,b,n,info)
           if (info.eq.0) then
c                                                
c             correction et test de la convergence
c
              residu = 0
              do 30 i=1,n
                x(i) = x(i) + b(i)
                residu = residu + b(i)*b(i)
30            continue  
              residu=sqrt(residu)       
              print*,' it = ',iter,' residu =',residu
              if (residu.lt.1.e-10) goto 50
           else
               print*,' pivot nul'
               stop
           end if
40      continue
c
c       resultats
c
50      continue       
        print*,' resultat'    
        do 60 i=1,n
             print*,' i =',i,' x(i) =',x(i)
60      continue
        end
c
c       la fonction c - c(x,0) - et sa derivee - c(x,1) - 
c
        double precision function c(x,n)
        implicit double precision (a-h,o-z)
        if (n.eq.0) then
                c=sin(x)
        else        
                c=cos(x)
        end if 
        end

Voilà le résultat du programme lorsque $\vec{x}_0=\vec{1}$ :

  iterations
  it =            1 residu =   4.599719439308855    
  it =            2 residu =   1.576327426641579    
  it =            3 residu =  0.1216083195869448    
  it =            4 residu =  8.0791194703510658E-05
  it =            5 residu =  2.9771670602180252E-14
  resultat
  i =           1 x(i) =  0.0000000000000000E+00
  i =           2 x(i) =  0.0000000000000000E+00
  i =           3 x(i) =  0.0000000000000000E+00
  i =           4 x(i) =  1.5777218104420236E-30
  i =           5 x(i) =  0.0000000000000000E+00
  i =           6 x(i) = -3.1554436208840472E-30
  i =           7 x(i) =  0.0000000000000000E+00
  i =           8 x(i) = -1.5777218104420236E-30
  i =           9 x(i) =  0.0000000000000000E+00
  i =          10 x(i) = -9.8607613152626476E-32