Versuchen Sie, das Programm "FORTRAN77 Numerical Computing Programming" auf C und Python zu portieren (Teil 2).

Vorheriger Artikel Portierung des Programms "FORTRAN77 Numerical Computing Programming" nach C und Python (Teil 1)

Aus dem Code von FORTRAN77 Numerical Calculation Programming, 1.2.

Anhäufung von Fehlern im Trapezgesetz

Kumulativer Rundungsfehler

Zitiert vom Ende von P.5 bis zur Vorderseite des P.6-Programms

Integration als Beispiel für einen kumulativen Rundungsfehler

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

 >, Während die Schrittweite $ h $ schrittweise verringert wird, ** Trapezregel **

>```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}

Versuchen Sie die numerische Integration mit>.

Fortran-Code

Ein Programm, das die numerische Integration von (1.14) und (1.15) durchführt.

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

Das Ausführungsergebnis ist

 ---- 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

Wird sein. Oh, die Ausgabe von ERR und L (ERR) unterscheidet sich vom Buch von ungefähr $ N = 128 $ ...

Den Code verstehen

Warum heißt dieses Programm überhaupt STRAP? Abgesehen von der Frage, was war die Integration? Lassen Sie uns einen zweiten Blick darauf werfen und überprüfen, was wir in das Buch und den Code schreiben. 　 In dem Buch wird die Integrationsgleichung durch das Trapezgesetz numerisch integriert. Viele Trapezoide werden in den Bereich des Graphen gestellt, der durch die ursprüngliche Integration dargestellt wird, und die Summe der Flächen der Trapezoide wird in den Integrationsgraphen integriert. Um es eine Annäherung an die Fläche von zu machen.

Das Programm führt also die Formeln (1.14) und (1.15) aus, während $ h $ immer kleiner wird. Das Thema ist, wie sehr es sich von der in (1.13) geforderten Antwort unterscheidet.

Ich habe ein Flussdiagramm gezeichnet.

Zunächst machen wir in * 1 $ arctan 1.0 $, das die inverse Dreiecksfunktion arctangent verwendet. Ausdrücken der Bogentangente als konstantes Integral,

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

Es wird die Formel. Wenn $ x $ $ 1 $ ist und $ z $ $ x $ ist, entspricht dies der Formel in (1.13), sodass EV hier der gleiche Wert wie (1.13) zugewiesen wird, dh $ I $. ..

Als nächstes wird in * 2 $ 2 $ durch $ 2 ^ K $ ersetzt, und $ N $ ist die Formel (1.14) für $ \ sum_ {j = 1} ^ {n-1} $. Die Anzahl der Bestellungen, die grob gesagt in dem Wert verwendet werden, der allmählich zunimmt. Das hier verwendete $ K $ ist ein Zähler für Schleife 1, und diese Schleife 1 dient dazu, das Fehlerverhalten bis zu 2 bis zur 13. Potenz zu testen.

Und für * 3 wurde dies ursprünglich als $ h = \ frac {ba} {n} $ in (1.14) ausgedrückt, und in (1.15) wurde es zu $$ a = 0, \ quad b = 1 . Wenn Sie es also ersetzen, ist es dasselbe wie $ h = \ frac {1} {n} $$.

Dann ist der erste Ausdruck in (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] Wenn Sie zuerst die zweite Formel anwenden $ 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]

Und wenn Sie (1.15) durch $ a = 0 $ und $ b = 1 $ ersetzen

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]

Und $ f (x) = \ frac {1} {1 + x ^ 2} $, also

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]

Von

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}

des Weiteren

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]

damit,

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}

Kann organisiert werden.

• 4 ersetzt $ S $ durch $ 0.5 \ times (1.0 + 0.5) , was $ \ frac {1} {2} \ times \ big (\ frac {2} {2} + ersetzt \ frac {1} {2} \ big) Wenn Sie darüber nachdenken, $ \ frac {1} {2} \ times \ frac {3} {2} = \ frac {3} {4} $ Und es ist der Wert, der in der Formel (r2) erscheint.

Der Zähler I von Schleife 2 ist $$ j , geschrieben in $ \ sum_ {j = 1} ^ {n-1} $$. Es wäre leicht zu verstehen, wenn Sie es J. geben würden.

Und in * 5 wird $ S + 1 \ div (1+ (H \ mal I) ^ 2) $ $ S $ zugewiesen, aber diese Formel lautet $ 1 \ div (1+ (H \) mal I) ^ 2) Wir werden den Teil von $ organisieren. Als allererstes $ 1 \ div1 + {\ big (\ frac {1} {n} \ times j \ big)} ^ 2 $ 　$\frac{1}{1+{\big(\frac{1}{n}\times j\big)}^2}$　　 Von $ \ frac {1} {1 + \ frac {j ^ 2} {n ^ 2}} $ und es erscheint in der Formel (r1).

Schleife 2 ist also (1.14)

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

Es wird derjenige sein, der den Teil von berechnet. In * 6 ist die erste Formel von (1.14) vervollständigt, und in * 7 die Differenz zur Formel von (1.13), dh die Differenz zwischen dem ursprünglich durch die Integralberechnung erhaltenen Berechnungsergebnis und dem Näherungswert unter Verwendung des Trapezgesetzes. Ich rechne. 　 Wenn in * 8 die in * 7 erhaltene Differenz ungleich 0 ist, wird der Absolutwert der Differenz berechnet ("ABS (ERR)") und der gemeinsame Logarithmus dieses Wertes berechnet ("ALOG10 ()"). .. Der gemeinsame Logarithmus ist der Wert von $, so dass $ x = 10 ^ a $ und in der Formel als $ a = \ log_ {10} x $$ ausgedrückt wird. Mit diesem Wert können Sie ein Diagramm wie folgt zeichnen: Diese Grafik zeigt, dass der kumulative Fehler proportional zu $ h ^ 2 $ abnimmt, aber nicht abnimmt, wenn sich der Fehler dem Computer Epsilon nähert.

Unterschied zu Büchern

Im Vergleich dazu hat sich der Wert von $ S $ von $ N = 128 $ leicht erhöht. Der Punkt, an dem die Abnahme des kumulativen Fehlers aufhört, ist etwas später. Der in diesem Buch verwendete Computer (darf nicht als Computer geschrieben werden) ist MV4000 AOS / VS, es handelt sich also um ein 16-Bit-Betriebssystem, nicht wahr? Es ähnelt also dem Diagramm, das ich ausgegeben habe. Das Diagramm wurde unter Verwendung der Gleitkommazahl mit doppelter Genauigkeit berechnet, die als Referenz im Buch angegeben ist. Es ist nur ähnlich, aber es ist ein bisschen anders, aber in der Buchumgebung ist Gleitkomma-Arithmetik eine 6-stellige Rundungsarithmetik. Ich denke, dies ist auf das Ergebnis einer Kürzungsoperation auf einem 64-Bit-Betriebssystem zurückzuführen (die Umgebung ist 32-Bit). Bitte lassen Sie mich wissen, wenn Sie ein Experte sind. Nachtrag: Ich habe mir in den Kommentaren wirklich einen Experten (@cure_honey) sagen lassen. Vielen Dank.

C-Code

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;
}

Das Ausführungsergebnis ändert sich nicht.

Python-Code

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)

#Zur Berechnung 1.0, 0.5, 0.0, -9.Gleitkommazahl mit einfacher Genauigkeit von 9
tmp = np.array([1.0, 0.5, 0.0, -9.9],dtype=np.float32)

#Liste zum Zeichnen von Grafiken
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))
    
    #Variablensatz für Grafik
    x.append(k+1)
    y.append(alerr)

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

Das Ausgabeergebnis ist das gleiche. Außerdem habe ich eine Grafik erstellt.

Zusammenfassung

Ich denke, dass der Python-Code sauberer sein wird. Und ich denke, es ist eine gute Idee, Numpy zu streuen, um es einfach zu machen. Möchten Sie die Genauigkeit des Fortran-Codes verdoppeln ... Nein, aber dann können die Ausführungsergebnisse nicht mit dem Buch verglichen werden. Muu.