Рассчитать два связанных уравнения с Numpy

Question

На этот раз это не очень важный вопрос, но, возможно, интересный.

Предположим, у нас есть две переменные x и y.Эти переменные зависят от времени (дискретное время).У нас есть начальные условия, и мы хотим их повторить со временем.Допустим, у нас есть x[0] = a и y[0] = b.Теперь мы хотим вычислить все заданные точки за небольшой промежуток времени и знаем следующее соотношение между этими двумя переменными:

x[n+1] = x[n] + y[n]
y[n+1] = y[n] + np.sin(x[n+1])

Конечно, мы можем сделать это с помощью цикла:

x[0], y[0] = a, b

for n in range(100): # just an arbitrary iteration
    x[n+1] = x[n] + y[n]
    y[n+1] = y[n] + np.sin(x[n+1])

Хорошо.Это возможно, если я не ошибся =).Что я хочу, так это, возможно, иметь гораздо лучший и более numpy -подобный способ решения проблемы без итерации.Я пытался придумать какие-то сдвиги или другие вещи.Я просто хочу вычислить массивы без цикла, потому что циклы действительно скучные.У меня только что была идея с рекурсивными вызовами функций, но я должен попробовать ее завтра утром.

Answer 1

Как сказали выше комментаторы, это не может быть красиво векторизовано. Если вы хотите быстрое исправление, которое должно помочь с скоростью, рассмотрите вариант inline из scipy.weave. Код ниже представляет собой скрипт, который дает рабочий пример этого для вашей желаемой дискретной системы. Это должно быть намного лучше, чем цикл for Python для больших массивов.

import scipy.weave as weave
import numpy as np

def simulate_x_and_y(x,y,a,b):

    x[0] = a; y[0] = b;
    n_range = len(x)

    code = """
    #include <math.h>

    for(int n = 1; n < n_range; n++){
        x(n) = x(n-1) + y(n-1);
        y(n) = y(n-1) + sin(x(n)); 
    }
    """

    weave.inline(code, ['x','y','n_range'], 
                 type_converters = weave.converters.blitz,
                 compiler = "gcc", headers=["<math.h>"]
                )


if __name__ == "__main__":

    x = np.zeros(10);
    y = np.zeros(10);
    a = 4.0; b = 10.0;

    print "Before:"
    print x
    print y

    simulate_x_and_y(x,y,a,b)

    print "After:"
    print x
    print y