Article précédent Portage du programme "FORTRAN77 Numerical Computing Programming" vers C et Python (Partie 1)

De Programmation de calcul numérique FORTRAN77, code 1.2.

Accumulation d'erreurs en loi trapézoïdale

Erreur d'arrondi cumulative

Cité de la fin de P.5 au début du programme P.6

Intégration comme exemple d'erreur d'arrondi cumulée

I = \int_{0}^{1}\frac{1}{1+x^2}dx = \frac{\pi}{4} = 0.7853982\tag{1.13}


 >, Tout en réduisant progressivement la largeur du pas $ h $, ** règle trapézoïdale **

>```math
I_h = h\biggr[\frac{1}{2}f(a)+\sum_{j=1}^{n-1}f(a+jh)+\frac{1}{2}f(b)\biggr], \quad　　h=\frac{b-a}{n}\tag{1.14}

a=0, \quad b=1, \quad f(x)=\frac{1}{1+x^2} \tag{1.15}

Essayez l'intégration numérique avec>.

Code Fortran

Un programme qui effectue l'intégration numérique de (1.14) et (1.15).

`strap.f`


*     Accumulation of round off error
      PROGRAM STRAP
        PARAMETER (ONE = 1.0)
*
        EV = ATAN(ONE)
*
        WRITE (*,2000)
2000    FORMAT (' ---- TRAP ----')
*
        DO 10  K = 1, 13
* 
            N = 2**K
            H = 1.0 / N
*
            S = 0.5 * (1.0 + 0.5)
            DO 20 I = 1, N - 1
                S = S + 1 / (1 + (H * I)**2)
  20        CONTINUE
            S = H * S
*
            ERR = S - EV
            IF (ERR .NE. 0.0) THEN
                ALERR = ALOG10 (ABS(ERR))
            ELSE
                ALERR = -9.9
            END IF
*
            WRITE (*,2001) N, H, S, ERR, ALERR
2001        FORMAT (' N=',I6,'  H=',1PE9.2,'  S=',E13.6,
     $              '  ERR=',E8.1,'  L(ERR)=',0PF4.1)
*
  10    CONTINUE
*
      END

Le résultat de l'exécution est

 ---- TRAP ----
 N=     2  H= 5.00E-01  S= 7.750000E-01  ERR=-1.0E-02  L(ERR)=-2.0
 N=     4  H= 2.50E-01  S= 7.827941E-01  ERR=-2.6E-03  L(ERR)=-2.6
 N=     8  H= 1.25E-01  S= 7.847471E-01  ERR=-6.5E-04  L(ERR)=-3.2
 N=    16  H= 6.25E-02  S= 7.852354E-01  ERR=-1.6E-04  L(ERR)=-3.8
 N=    32  H= 3.12E-02  S= 7.853574E-01  ERR=-4.1E-05  L(ERR)=-4.4
 N=    64  H= 1.56E-02  S= 7.853879E-01  ERR=-1.0E-05  L(ERR)=-5.0
 N=   128  H= 7.81E-03  S= 7.853957E-01  ERR=-2.4E-06  L(ERR)=-5.6
 N=   256  H= 3.91E-03  S= 7.853976E-01  ERR=-5.4E-07  L(ERR)=-6.3
 N=   512  H= 1.95E-03  S= 7.853978E-01  ERR=-4.2E-07  L(ERR)=-6.4
 N=  1024  H= 9.77E-04  S= 7.853981E-01  ERR=-6.0E-08  L(ERR)=-7.2
 N=  2048  H= 4.88E-04  S= 7.853982E-01  ERR= 6.0E-08  L(ERR)=-7.2
 N=  4096  H= 2.44E-04  S= 7.853983E-01  ERR= 1.2E-07  L(ERR)=-6.9
 N=  8192  H= 1.22E-04  S= 7.853981E-01  ERR=-6.0E-08  L(ERR)=-7.2

Sera. Oh, la sortie de ERR et L (ERR) est différente du livre d'environ $ N = 128 $ ...

Comprendre le code

Pourquoi le nom de ce programme STRAP en premier lieu? Mis à part la question, quelle a été l'intégration? Jetons un second regard et vérifions ce que nous écrivons dans le livre et le code. 　 Dans le livre, l'équation d'intégration est numériquement intégrée par la loi trapézoïdale. De nombreux trapèzes sont placés dans la zone du graphe représentée par l'intégration d'origine, et la somme des aires des trapèzes est intégrée dans le graphe d'intégration. Pour en faire une approximation de la superficie de.

Donc, ce que fait le programme est d'exécuter les formules (1.14) et (1.15) tout en rendant $ h $ de plus en plus petit. Le thème est combien il diffère de la réponse requise dans (1.13).

J'ai dessiné un organigramme.

Tout d'abord, ce que nous faisons dans * 1 est $ arctan 1.0 $, qui utilise la fonction triangulaire inverse arctangente. Exprimer l'arc tangent comme une intégrale constante,

arctan\,x = \int_0^x \frac 1 {z^2 + 1}\,dz

Cela devient la formule. Si $ x $ est $ 1 $ et $ z $ est $ x $, ce sera la même chose que la formule dans (1.13), donc la même valeur que (1.13), c'est-à-dire que $ I $ sera affecté à EV ici. ..

Ensuite, dans * 2, $ 2 ^ K $ est remplacé par $ N $, et $ N $ est la formule (1.14) de $ \ sum_ {j = 1} ^ {n-1} $. Le nombre de commandes utilisées dans, en gros, la valeur qui augmente progressivement. Le $ K $ utilisé ici est un compteur de boucle 1, et cette boucle 1 consiste à essayer le comportement d'erreur jusqu'à 2 à la 13e puissance.

Et pour * 3, cela a été initialement exprimé par $ h = \ frac {ba} {n} $ dans (1.14), et dans (1.15) il est devenu $$ a = 0, \ quad b = 1 . Donc, si vous le remplacez, c'est la même chose que $ h = \ frac {1} {n} $$.

Ensuite, la première expression dans (1.14) $ I_h = h \ biggr [\ frac {1} {2} f (a) + \ sum_ {j = 1} ^ {n-1} f (a + jh) ) + \ Frac {1} {2} f (b) \ biggr] Si vous appliquez d'abord la deuxième formule $ h = \ frac {1} {n} $$

I_\frac{1}{n}=\frac{1}{n}\biggr[\frac{1}{2}f(a)+\sum_{j=1}^{n-1}f(a+j\frac{1}{n})+\frac{1}{2}f(b)\biggr]

Et si vous remplacez (1.15) par $ a = 0 $ et $ b = 1 $

I_\frac{1}{n}=\frac{1}{n}\biggr[\frac{1}{2}f(0)+\sum_{j=1}^{n-1}f(0+\frac{j}{n})+\frac{1}{2}f(1)\biggr]

Et $ f (x) = \ frac {1} {1 + x ^ 2} $, donc

I_\frac{1}{n}=\frac{1}{n}\biggr[\frac{1}{2}\times\frac{1}{1+0^2}+\sum_{j=1}^{n-1}\frac{1}{1+\big(\frac{j}{n}\big)^2} +\frac{1}{2}\times\frac{1}{1+1^2}\biggr]

I_\frac{1}{n}=\frac{1}{n}\biggr[\frac{1}{2}+\sum_{j=1}^{n-1}\frac{n}{1+\frac{j^2}{n^2}}+\frac{1}{4}\biggr]\tag{r1}

plus loin

I_\frac{1}{n}=\frac{1}{n}\biggr[\frac{1}{2}+\sum_{j=1}^{n-1}\frac{n^2}{n^2+j^2}+\frac{1}{4}\biggr]

alors,

I_\frac{1}{n}=\frac{1}{n}\biggr[\frac{3}{4}+\sum_{j=1}^{n-1}\frac{n^2}{n^2+j^2}\biggr]\tag{r2}

Peut être organisé.

4 substitue $ 0.5 \ times (1.0 + 0.5) $ pour $ S , qui est remplacé par $ \ frac {1} {2} \ times \ big (\ frac {2} {2} + \ frac {1} {2} \ big) Si vous y réfléchissez, $ \ frac {1} {2} \ times \ frac {3} {2} = \ frac {3} {4} $ Et c'est la valeur qui apparaît dans la formule (r2).

Le compteur I de la boucle 2 est $$ j $ écrit en $ \ sum_ {j = 1} ^ {n-1} $. Ce serait facile à comprendre si vous le donniez à J.

Et dans * 5, $ S + 1 \ div (1+ (H \ times I) ^ 2) $ est assigné à $ S , mais cette formule est $ 1 \ div (1+ (H ) fois I) ^ 2) Nous organiserons la partie de $. Tout d'abord $ 1 \ div1 + {\ big (\ frac {1} {n} \ times j \ big)} ^ 2 $ 　$\frac{1}{1+{\big(\frac{1}{n}\times j\big)}^2}$　　 De $ \ frac {1} {1 + \ frac {j ^ 2} {n ^ 2}} $$ et il apparaît dans la formule (r1).

Donc, la boucle 2 est (1.14)

\frac{1}{2}f(a)+\sum_{j=1}^{n-1}f(a+jh)+\frac{1}{2}f(b)

Ce sera celui qui calcule la part de. En * 6, la première équation de (1.14) est complétée, et en * 7, la différence avec l'équation de (1.13), c'est-à-dire la différence entre le résultat du calcul obtenu à l'origine par le calcul intégral et la valeur approximative utilisant la loi trapézoïdale. Je calcule. 　 En * 8, lorsque la différence obtenue en * 7 n'est pas égale à 0, la valeur absolue de la différence est calculée (ʻABS (ERR) ), et la valeur logarithmique commune de cette valeur est calculée (ʻALOG10 ()). .. Le logarithme courant est la valeur de $ telle que $ x = 10 ^ a $, qui est exprimée par $ a = \ log_ {10} x $$ dans la formule. Vous pouvez utiliser cette valeur pour tracer un graphique comme celui-ci: Ce graphique montre que l'erreur cumulée diminue proportionnellement à $ h ^ 2 $, mais elle ne diminue pas à mesure que l'erreur s'approche de l'ordinateur d'Epsilon.

Différence avec les livres

Par comparaison, la valeur de $ S $ a légèrement augmenté de $ N = 128 $. Le point où la diminution de l'erreur cumulée s'arrête est un peu plus tard. L'ordinateur utilisé dans le livre (n'ose pas être écrit comme ordinateur) est MV4000 AOS / VS, donc c'est un système d'exploitation 16 bits, n'est-ce pas? Donc, c'est similaire au graphique que j'ai publié, le graphique calculé en utilisant le nombre à virgule flottante double précision qui est donné comme référence dans le livre. C'est juste similaire, mais c'est un peu différent, mais dans l'environnement du livre, l'arithmétique en virgule flottante est une arithmétique d'arrondi à 6 chiffres. Donc, je suppose que cela est dû au résultat d'une opération de troncature sur un système d'exploitation 64 bits (l'environnement est 32 bits). Veuillez me faire savoir si vous êtes un expert. Postscript: J'ai vraiment eu un expert (@cure_honey) me le dire dans les commentaires. Merci beaucoup.

Code C

`strap.c`


#include <stdio.h>
#include <math.h>

int main(void)
{
    const float ONE = 1.0f;
    float ev = atanf(ONE);
    int k, i;
    int n;
    float h, s, err, alerr;

    printf("%14.7E\n", ev);

    printf("--- TRAP ---\n");

    for(k=1; k<=13; k++)
    {
        n = pow(2,k);
        h = 1.0f / n;
        s = 0.5f * (1.0f + 0.5f);
        
        for(i=1; i<=n-1; i++)
        {
            s = s + 1 / (1+ powf(h*i,2));
        }
        
        s = h * s;
        err = s - ev;

        if(err != 0.0f)
        {
            alerr = log10(fabsf(err));
        }
        else 
        {
            alerr = -9.9f;
        }
        printf("N=%6d　H=%9.2E　S=%13.6E　ERR=%8.1E　L(ERR)=%4.1F\n",n, h, s, err, alerr);
    }

    return 0;
}

Le résultat de l'exécution ne change pas.

Code Python

`strap.py`


import numpy as np
import matplotlib.pyplot as plt

ONE = np.array((1.0),dtype=np.float32)
ev = np.arctan(ONE,dtype=np.float32)

h, s, err, alerr = np.array([0.0, 0.0, 0.0, 0.0],dtype=np.float32)

#Pour le calcul 1.0, 0.5, 0.0, -9.Nombre à virgule flottante simple précision de 9
tmp = np.array([1.0, 0.5, 0.0, -9.9],dtype=np.float32)

#Liste pour dessiner des graphiques
x = []
y = []

print("--- TRAP ---")

for k in range(13):
    n = np.power(2,k+1)
    h = tmp[0] / n.astype(np.float32)
    s = tmp[1] * (tmp[0] + tmp[1])
    
    for i in range(n-1):
        s = s + tmp[0] / (tmp[0] + np.square(h*(i+1),dtype=np.float32))  
        
    s = s * h
    err = s - ev

    if err != tmp[2]:
        alerr = np.log10(np.abs(err))
    else:
        alerr = tmp[3]

    print("N=%6d　H=%9.2E　S=%13.6E　ERR=%8.1E　L(ERR)=%4.1F" % (n, h, s, err, alerr))
    
    #Jeu de variables pour le graphique
    x.append(k+1)
    y.append(alerr)

#Dessin graphique
fig, ax = plt.subplots()
ax.scatter(x,y)
ax.set_xlabel("n (1/2^n)")
ax.set_ylabel("l (10^l)")
plt.show()

Le résultat de sortie est le même. Aussi, j'ai mis un graphique.

Résumé

Je pense que le code Python va être plus propre. Et je pense que c'est une bonne idée de disperser numpy pour faire simple. Voulez-vous doubler la précision du code Fortran ... Non, mais alors les résultats d'exécution ne peuvent pas être comparés avec le livre. Muu.

Essayez de porter le programme "FORTRAN77 Numerical Computing Programming" vers C et Python (partie 2)