Почему моя реализация метода Рунге-Кутты четвертого порядка не работает?

После поиска на этом сайте и в моем справочнике я узнал, что понятия не имею, почему мой код не работает.Почему моя реализация метода Рунге-Кутты четвертого порядка не работает?

Я сделал реализацию Рунге-Кутты четвертого порядка для системы массовых весов (с амортизацией), так как мой профессор показал нам класс. Однако, как видно, получившийся граф довольно странный. example image

код я в конечном итоге написание это:

#! /usr/bin/env python3 
#-*- coding: utf-8 -*- 

def f(data, t, x1, v1): 
    from math import cos 

    F = data["F"] 
    c = data["c"] 
    k = data["k"] 
    m = data["m"] 
    omega = data["omega"] 

    return([v1, (F*cos(omega*t) - c*v1 - k*x1)/m]) 

def run(data = {}): 
    xi, vi, ti = [data["u1"]], [data["v1"]], [data["t_ini"]] 
    dt = data["dt"] 

    while ti[-1] <= data["t_end"]: 
     xn = xi[-1] 
     vn = vi[-1] 
     tn = ti[-1] 

     K1 = f(data, t = tn, x1 = xn, v1 = vn) 
     K1 = [dt*K1[i] for i in range(len(K1))] 

     K2 = f(data, t = tn + 0.5*dt, x1 = xn + 0.5*K1[0], v1 = vn + 0.5*K1[1]) 
     K2 = [dt*K2[i] for i in range(len(K2))] 

     K3 = f(data, t = tn + 0.5*dt, x1 = xn + 0.5*K2[0], v1 = vn + 0.5*K2[1]) 
     K3 = [dt*K3[i] for i in range(len(K3))] 

     K4 = f(data, t = tn + dt, x1 = xn + K3[0], v1 = vn + K3[1]) 
     K4 = [dt*K4[i] for i in range(len(K4))] 

     xn = xn + (K1[0] + 2*K2[0] + 2*K3[0] + K4[0])/6 
     vn = xn + (K1[1] + 2*K2[1] + 2*K3[1] + K4[1])/6 

     ti.append(tn+dt) 
     xi.append(xn) 
     vi.append(vn) 

    return(ti, xi, vi)

Это импортируется файл main.py, который содержит только с графическим интерфейсом и сюжетные части программы, а функция была выведена в классе , поэтому я считаю, что ошибка находится в самой Рунге-Кутте. (Наверное, это что-то глупое, что я испортил.)

Я попытался переключить K в «xn» и «vn», грубые форсирующие значения «F» и «c» в f(), переписывая все и записывая каждый элемент каждого K вручную (как K11, K12, K21 и т. д.), но он дает только экспоненциальный результат. Кроме того, включение возврата f() для массива numpy ничего не помогло.

Я нашел несколько вопросов для метода RK4, но я не мог решить эту проблему и не понял, что случилось. У меня есть некоторый knoledge метода, но это на самом деле моя первая реализация, поэтому, пожалуйста, любая помощь приветствуется.

Я использую дистрибутив Anaconda для python3, если это имеет значение.

источник

2013-11-30 Ryu

http://scicomp.stackexchange.com/? –

Я думал, что это больше связано с python, чем с самим rk ... например, испортить какой-то список или передать неправильный аргумент для функции – Ryu

ok np (похоже, что ответы, которые здесь отвечают, меняются, с меньшими ответами на этот вид вопрос, поэтому я подумал, что это может помочь). –

Может ли быть эта линия?

vn = xn + (K1[1] + 2*K2[1] + 2*K3[1] + K4[1])/6

Должна ли она быть:

vn = vn + (K1[1] + 2*K2[1] + 2*K3[1] + K4[1])/6

источник

2013-11-30 22:47:05 dansalmo

Человек, я чувствую себя таким глупым сейчас ... Большое вам спасибо! – Ryu

Чувство дурака - это регулярная часть моего опыта программирования. – dansalmo

def rk4(f, xvinit, Tmax, N): 
T = np.linspace(0, Tmax, N+1) 
xv = np.zeros((len(T), len(xvinit))) 
xv[0] = xvinit 
h = Tmax/N 
for i in range(N): 
    k1 = f(xv[i]) 
    k2 = f(xv[i] + h/2.0*k1) 
    k3 = f(xv[i] + h/2.0*k2) 
    k4 = f(xv[i] + h*k3) 
    xv[i+1] = xv[i] + h/6.0 *(k1 + 2*k2 + 2*k3 + k4) 
return T, xv

Это то, что вы можете использовать для выполнения функции с помощью алгоритма RK4

источник

2014-12-07 16:46:27

Почему моя реализация метода Рунге-Кутты четвертого порядка не работает?

ответ

Смежные вопросы