[Wissenschaftlich-technische Berechnung nach Python] Numerische Integration, Trapezgesetz / Simpson-Gesetz, numerische Berechnung, scipy

Die numerische Integration diskreter Daten erfolgt mit der Cumtrapz-Methode (Trapezgesetz) und der Simps-Methode (Simpson-Gesetz) von scipy.integrate. Betrachten Sie als Beispiel $ \ int_0 ^ 1 \ frac {4} {1 + x ^ 2} dx = \ pi $.

Inhalt

(1.A) Code mit scipy. ** Dies ist, wenn Sie es eilig haben. ** ** ** (2) Nachtrag zur Berechnungsgenauigkeit des Trapezgesetzes und des Simpsonschen Gesetzes. (3) Einfache Implementierung von Trapez- und Simpson-Regeln mit Python-Code. Nützlich, wenn Sie das Berechnungsverfahren kennenlernen möchten.

(1) Code mit scipy

(1) Scipy-Verwendungscode. Prägnant.

from scipy import integrate
import numpy as np

x=np.linspace(0,1,5)  # [0,1]Speichern Sie das in 5 gleiche Teile unterteilte Gitter in x
y=4/(1+x**2)          #Integrand der numerischen Integration

y_integrate_trape = integrate.cumtrapz(y, x)  #Numerische Integralberechnung nach dem Trapezgesetz
y_integrate_simps = integrate.simps(y, x) #Numerische Integralberechnung nach dem Simpsonschen Gesetz

print(y_integrate_trape[-1]) #Ergebnisse anzeigen
print(y_integrate_simps)     #Ergebnisse anzeigen

Ergebnis (1): Scipy verwenden

3.13117647059   #Trapezregel
3.14156862745   #Simpsons Gesetz
3.14159265358979...Genauer Wert(π)

(2) [Nachtrag] Vergleich der Genauigkeit und Konvergenzgeschwindigkeit des Trapezgesetzes und des Simpsonschen Gesetzes

Bekanntlich ist für eine ausreichend kleine Schrittweite h, Der Integrationsfehler nach dem Trapezgesetz beträgt $ O (h ^ 2) $, und der des Simpson-Gesetzes beträgt $ O (h ^ 3) $. Unter der Annahme, dass die Anzahl der Gitter N ist, sind dies $ O (N ^ {-2}) $ bzw. $ O (N ^ {-3}) $. Es ist lehrreich, dies direkt anhand dieses Beispiels zu überprüfen.

Unten finden Sie den Code zur Bestätigung. Beide logarithmischen Diagramme werden mit dem relativen Fehler aufgrund der numerischen Integration auf der y-Achse und der Anzahl der Gitter N auf der horizontalen Achse erstellt. In dem Bereich, in dem die obige Beziehung gilt, entspricht $ log y \ propto n logN $, $ n = -2 $ dem Trapezgesetz und $ n = -3 $ dem Simpson-Gesetz.

(2)Vergleich der Berechnungsgenauigkeit

from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt

err_trape=[]
err_simps=[]
nn=[4,8,16,32,64,128,256,512,1024,2048] #Speichern Sie die Anzahl der Raster von 4 bis 2048 in Liste nn
for j in nn:
    x=np.linspace(0,1, j)
    y=4/(1+x**2)
    y_integrate_trape = integrate.cumtrapz(y, x) #Numerische Integralberechnung nach dem Trapezgesetz
    y_integrate_simps = integrate.simps(y, x)     #Numerische Integralberechnung nach dem Simpsonschen Gesetz
    err_trape.append(abs(np.pi-y_integrate_trape[-1])/np.pi)  #Relative Fehlerbewertung der numerischen Integration nach dem Trapezgesetz
    err_simps.append(abs(np.pi-y_integrate_simps)/np.pi)     #Relative Fehlerbewertung der numerischen Integration nach dem Simpsonschen Gesetz

# for plot
ax = plt.gca()
ax.set_yscale('log')  #Zeichnen Sie die y-Achse auf der Protokollskala
ax.set_xscale('log')  #Zeichnen Sie die x-Achse auf der Protokollskala
plt.plot(nn,err_trape,"-", color='blue', label='Trapezoid rule')
plt.plot(nn,err_simps,"-", color='red', label='Simpson rule')


plt.xlabel('Number of grids',fontsize=18)
plt.ylabel('Error (%)',fontsize=18)
plt.grid(which="both") #Rasteranzeige."both"Zeichnet ein Gitter auf beiden xy-Achsen.

plt.legend(loc='upper right')


plt.show()

Die Steigung der geraden Linie in der Figur entspricht $ n $ der Konvergenzreihenfolge. Wie erwartet ist das Trapezgesetz (blaue Linie) $ n \ sim-2 $ und das Simpson-Gesetz (rote Linie) ist $ n \ sim-3 $.

(3) Einfache Implementierung von Trapez- und Simpson-Regeln mit Python-Code.

Verwendung: Bereiten Sie eine XY-Datendatei vor ('xy.dat')

python3 fuga.py xy.dat

Anschließend wird das numerische Berechnungsergebnis nach der Trapezregel und der Simpson-Regel angezeigt.

fuga.py

import os, sys, math


def integrate_trape(f_inp): #Funktion der numerischen Integration nach dem Trapezgesetz
    x_lis=[]; y_lis=[]
#   f_inp.readline()
    for line in f_inp:
        x_lis.append(float(line.split()[0]))
        y_lis.append(float(line.split()[1]))
    sum_part_ylis=sum(y_lis[1:-2]) 
    del_x=(x_lis[1]-x_lis[0])
    integrated_value = 0.5*del_x*(y_lis[0]+y_lis[-1]+2*sum_part_ylis) 
    print("solution_use_trapezoid_rule: ", integrated_value)

def integrate_simpson(f_inp): #Funktion der numerischen Integration nach dem Simpsonschen Gesetz
    x_lis=[]; y_lis=[]
    n_counter = 0
    for line in f_inp:
        x_lis.append(float(line.split()[0]))
        if n_counter % 2 == 0:
            y_lis.append(2*float(line.split()[1]))
        else:
            y_lis.append(4*float(line.split()[1]))
        n_counter += 1
    sum_part_ylis=sum(y_lis[1:-2]) # sum from y_1 to y_n-1 
    del_x=(x_lis[1]-x_lis[0])
    integrated_value = del_x*(y_lis[0]+y_lis[-1]+sum_part_ylis)/3 
    print("solution_use_Simpson_formula: ", integrated_value)

##
def main():
    fname=sys.argv[1]

    f_inp=open(fname,"r")
    integrate_trape(f_inp) 
    f_inp.close()

    f_inp=open(fname,"r")
    integrate_simpson(f_inp) 
    f_inp.close()

if __name__ == '__main__':
    main()