Rozwiązuję numerycznie dla x (t) dla układu równań różniczkowych pierwszego rzędu. System jest:Metoda Vectorize Forward Euler dla układu równań różniczkowych
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
Wersja "wektoryzowana" jedynie wektoryzuje "h * f (x [step])" lub tylko dwie operacje. Dodatkowy koszt utworzenia tablicy numpy kompensuje wzrost prędkości. W zależności od tego, co robisz, możesz zajrzeć do [scipy.integrate.ode] (http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html#scipy.integrate. oda). – Daniel