Рунге – Кутта 4-го порядка

Давайте рассмотрим, что у меня есть система из 4 ОДУ: dX / dt = F (X), где X - вектор (4-размерность), а F: R ^ 4 -> R ^ 4. F называется vectorDE_total_function, и я пытаюсь вычислить решение с помощью RK-4.

def solvingDES():
    previous_vector = np.array ([theta_1, omega_1, theta_2, omega_2]);
    for current_time in time:
        temp_vector = previous_vector;
        RK_vector = np.array([0.0,0.0,0.0,0.0]);
        for c in [6,3,3,6]:
            RK_vector = vectorDE_total_function(previous_vector + c * RK_vector/6) * time_step;
            temp_vector += RK_vector / c;
        previous_vector = temp_vector;
        current_time += 1;

И похоже, что где-то я ошибаюсь, но не знаю где. Это кажется законным?


python numpy numerical-methods runge-kutta
person Mikhail Tikhonov    schedule 11.12.2016    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Это странный, редко встречающийся способ его реализации, и он работает только для классического RK4, другие методы Рунге-Кутты не будут работать так. Но общая идея кажется верной.

У вас обычная ошибка в обычно неожиданном месте. Параметр

temp_vector = previous_vector;

и позже

previous_vector = temp_vector;

вы не делаете копии векторов, но заставляете обе объектные ссылки использовать одни и те же векторы. Использовать

temp_vector = previous_vector.copy();

or

previous_vector = temp_vector[:];

для принудительного копирования векторных данных.

person Lutz Lehmann    schedule 11.12.2016
comment
Тогда я не уверен, что делает '='. Является ли это обязательным, а не поручением? - person Mikhail Tikhonov; 11.12.2016
comment
Векторные переменные являются ссылками или указателями на фактические объекты памяти. Таким образом, = копирует указатели. Возможно, но не обязательно, с некоторым подсчетом ссылок. - += также воздействует на объект, на который указывает на левой стороне, таким образом, в отличие от a=a+b, где новый a является новым объектом, в a+=b a это тот же объект, что и раньше. Таким образом, temp_vector += RK_vector / c; не нарушает совпадение указателей и, таким образом, также изменяет previous_vector (… потому что никаких дальнейших действий не происходит). - person Lutz Lehmann; 11.12.2016
comment
Так что нет, это не привязка, поскольку она поддерживает синхронизацию двух объектов с помощью некоторого механизма обработки событий, то есть изменения в объекте A, вызывающего вызов метода объекта B. Ничего подобного здесь нет. - person Lutz Lehmann; 11.12.2016
comment
Я в основном использую C-код, так что .. Вы утверждаете, что что-то вроде b = a, a + = b также вызывает изменение b? Я думал, что b передается по ссылке (т.е. не точный адрес памяти, а значение) - person Mikhail Tikhonov; 11.12.2016
comment
Да, у меня именно это произошло в каком-то коде RK4. ylast=ycurr; ycurr+=h*k;, где k - усредненный наклон шага RK, а затем задался вопросом, почему в результирующем графике были эти странные колебания при интерполяции от последнего до текущего значения. Для правильного поведения требовалось явное добавление ycurr=ycurr+h*k;. - person Lutz Lehmann; 11.12.2016
comment
Код RK4 с адаптивным временным шагом и кубической интерполяцией Безье находится здесь: stackoverflow.com/a/40866829/3088138 - person Lutz Lehmann; 11.12.2016
comment
Немного помогло, вот осцилляции появились, но все равно странные. Я предполагаю, что это код проблемы с F-функцией. Спасибо, попробую исправить. - person Mikhail Tikhonov; 11.12.2016