Metoda Vectorize Forward Euler dla układu równań różniczkowych

Question

Rozwiązuję numerycznie dla x (t) dla układu równań różniczkowych pierwszego rzędu. System jest:

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

I wprowadziły metodę Eulera przodu rozwiązać ten problem w następujący sposób:

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 

teraz chciałbym wektoryzacji kod dalej i utrzymać x i y w taki sam tablicy i pochodzą z następującego roztworu:

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 

The pytanie: forward_euler_vector() działa, ale czy to było najlepsze, aby go wektoryzować? Pytam, ponieważ wersja wektorowy biegnie około 20 ms wolniejszy na moim laptopie:

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 komentarz wyjaśnia bardzo dobrze, co się dzieje. Wywołanie array() w ramach f(x) wprowadza pewne obciążenie, które zabija korzyść użycia mnożenia macierzy w wyrażeniu h * f(x[step]).

A jak mówi, być może zainteresuje Cię zestaw scipy.integrate dla ładnego zestawu integratorów numerycznych.

Aby rozwiązać problem związany z wektorowaniem kodu, nie należy odtwarzać tablicy za każdym razem, gdy wywoływana jest f. Chciałbyś zainicjować tablicę raz i powrócić ją zmodyfikowaną przy każdym wywołaniu. Jest to podobne do zmiennej static w C/C++.

Można to osiągnąć za pomocą zmiennego argumentu domyślnego, który jest interpretowany tylko raz, w momencie definicji funkcji f(x) i ma zakres lokalny. Ponieważ ma być zmienny, należy ująć ją w liście jednego elementu:

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] 

Dzięki tej modyfikacji do kodu, napowietrznej tworzenia tablicy znika, a na moim komputerze mam zdobyć niewielką poprawę:

%timeit forward_euler()  #258ms 
%timeit forward_euler_vector() #248ms 

Oznacza to, że przyrost optymalizacji mnożenia macierzy za pomocą numpy jest niewielki, przynajmniej pod względem problemu.

Możesz od razu odłączyć się od funkcji f, wykonując operacje w pętli for, pozbawiając się narzutów związanych z połączeniem. Ta sztuczka domyślnego argumentu może być jednak zastosowana także z bardziej globalnymi integratorami czasu, gdzie musisz podać funkcję f.

EDIT: jak wskazał Jaime, innym sposobem jest w leczeniu static_tmp jako atrybut funkcji f i stworzyć go po zgłoszeniu funkcja ale przed wywołaniem go:

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

Nie zawsze jest banalna autojit rozwiązanie:

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