Vectorize Forward Euler-Methode für das System der Differentialgleichungen

Question

Ich bin numerisch für x (t) für ein System der Differentialgleichungen erster Ordnung zu lösen. Das System ist:

dx/dt = y
dy/dt = x - a * y (x^2 + y^2 -1)

I der Vorwärts-Euler-Verfahren implementiert haben, dies zu lösen, Problem wie folgt:

def forward_euler(): 
    h = 0.01 
    num_steps = 10000 

    x = np.zeros([num_steps + 1, 2]) # steps, number of solutions 
    y = np.zeros([num_steps + 1, 2]) 
    a = 1. 

    x[0, 0] = 10. # initial condition 1st solution 
    y[0, 0] = 5. 

    x[0, 1] = 0. # initial condition 2nd solution 
    y[0, 1] = 0.0000000001 

    for step in xrange(num_steps): 
     x[step + 1] = x[step] + h * y[step] 
     y[step + 1] = y[step] + h * (-x[step] - a * y[step] * (x[step] ** 2 + y[step] ** 2 - 1)) 

    return x, y 

Nun möchte ich mag den Code weiter vektorisieren und im gleichen Array x und y zu halten, habe ich mit der folgenden Lösung kommen:

def forward_euler_vector(): 
    num_steps = 10000 
    h = 0.01 

    x = np.zeros([num_steps + 1, 2, 2]) # steps, variables, number of solutions 
    a = 1. 

    x[0, 0, 0] = 10. # initial conditions 1st solution 
    x[0, 1, 0] = 5. 

    x[0, 0, 1] = 0. # initial conditions 2nd solution 
    x[0, 1, 1] = 0.0000000001 

    def f(x): 
     return np.array([x[1], 
         -x[0] - a * x[1] * (x[0] ** 2 + x[1] ** 2 - 1)]) 

    for step in xrange(num_steps): 
     x[step + 1] = x[step] + h * f(x[step]) 

    return x 

des Frage: forward_euler_vector() funktioniert, aber war das der beste Weg, um es zu vektorisieren? Ich frage, weil die vektorisierte Version läuft etwa 20 ms langsamer auf meinem Laptop:

In [27]: %timeit forward_euler() 
1 loops, best of 3: 301 ms per loop 

In [65]: %timeit forward_euler_vector() 
1 loops, best of 3: 320 ms per loop 

Answer 1

@Ophion Kommentar erklärt sehr gut, was los ist. Der Aufruf an array() innerhalb von f(x) führt zu etwas Overhead, was den Vorteil der Verwendung der Matrixmultiplikation im Ausdruck h * f(x[step]) zunichte macht.

Und wie er sagt, Sie können interessiert sein, einen Blick auf scipy.integrate für eine schöne Reihe von numerischen Integratoren.

Um das Problem bei der Vektorisierung Ihres Codes zu lösen, möchten Sie das Array nicht jedes Mal neu erstellen, wenn Sie f anrufen. Sie möchten das Array einmalig initialisieren und bei jedem Aufruf modifizieren. Dies ist vergleichbar mit einer Variablen static in C/C++.

Sie können dies mit einem veränderbaren Standardargument erreichen, das zum Zeitpunkt der Definition der Funktion f(x) einmal interpretiert wird und einen lokalen Gültigkeitsbereich hat.

def f(x,static_tmp=[empty((2,2))]): 
    static_tmp[0][0]=x[1] 
    static_tmp[0][1]=-x[0] - a * x[1] * (x[0] ** 2 + x[1] ** 2 - 1) 
    return static_tmp[0] 

Mit dieser Änderung an Ihrem Code, der Aufwand für Array-Erstellung verschwindet, und auf meinem Rechner erhalte ich eine kleine Verbesserung: Da es wandelbar sein muss, können Sie es in einer Liste eines einzelnen Elements kapseln

Dies bedeutet, dass die Optimierung der Matrixmultiplikation mit numpy relativ gering ist, zumindest in Bezug auf das vorliegende Problem.

Sie können auch sofort die Funktion f loswerden, ihre Operationen innerhalb der for-Schleife tun, den Anruf Overhead loswerden. Dieser Trick des Standardarguments kann jedoch auch mit scipy allgemeineren Zeitintegratoren angewendet werden, wo Sie eine Funktion f bereitstellen müssen.

EDIT: wie von Jaime wies darauf hin, einen anderen Weg zu gehen static_tmp als Attribut der Funktion f, zu behandeln und zu erstellen, nachdem die Funktion deklariert haben, aber bevor sie rufen:

def f(x): 
    f.static_tmp[0]=x[1] 
    f.static_tmp[1]=-x[0] - a * x[1] * (x[0] ** 2 + x[1] ** 2 - 1) 
    return f.static_tmp 
f.static_tmp=empty((2,2)) 

Answer 2

Es ist immer die triviale Lösung autojit:

def forward_euler(initial_x, initial_y, num_steps, h): 

    x = np.zeros([num_steps + 1, 2]) # steps, number of solutions 
    y = np.zeros([num_steps + 1, 2]) 
    a = 1. 

    x[0, 0] = initial_x[0] # initial condition 1st solution 
    y[0, 0] = initial_y[0] 

    x[0, 1] = initial_x[1] # initial condition 2nd solution 
    y[0, 1] = initial_y[1] 

    for step in xrange(int(num_steps)): 
     x[step + 1] = x[step] + h * y[step] 
     y[step + 1] = y[step] + h * (-x[step] - a * y[step] * (x[step] ** 2 + y[step] ** 2 - 1)) 

    return x, y 

Timings:

from numba import autojit 
jit_forward_euler = autojit(forward_euler) 

%timeit forward_euler([10,0], [5,0.0000000001], 1E4, 0.01) 
1 loops, best of 3: 385 ms per loop 

%timeit jit_forward_euler([10,0], [5,0.0000000001], 1E4, 0.01) 
100 loops, best of 3: 3.51 ms per loop