Проблема с методом Рунге-Кутты для закона Кулона - PullRequest
Новая
       45

Проблема с методом Рунге-Кутты для закона Кулона

0 голосов
/ 22 мая 2019

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

Начиная с положительной начальной скорости vx0 и отрицательной начальной x-позиции x0, код работает нормально. Но когда я переворачиваю их так, чтобы vx0 было отрицательным, а x0 - положительным, я получаю разные ответы, когда они должны быть симметричными.

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

enter image description here

Вот где я вычисляю константы: 

# Initial conditions
time[0] = t = t0
vx1[0] = vx = vx0
vy1[0] = vy = vy0
x1[0] = x = x0
y1[0] = y = y0

for i in range(1, n + 1):
    # Calculate our constants
    k1vx = step * ax(x, y, q1, q2, m)
    k1vy = step * ay(x, y, q1, q2, m)
    k1x = step * vx
    k1y = step * vy
    k2vx = step * ax(x + 0.5 * k1x, y + 0.5 * k1y, q1, q2, m)
    k2vy = step * ay(x + 0.5 * k1x, y + 0.5 * k1y, q1, q2, m)
    k2x = step * (vx + (k1vx / 2))
    k2y = step * (vy + (k1vy / 2))
    k3vx = step * ax(x + 0.5 * k2x, y + 0.5 * k2y, q1, q2, m)
    k3vy = step * ay(x + 0.5 * k2x, y + 0.5 * k2y, q1, q2, m)
    k3x = step * (vx + (k2vx / 2))
    k3y = step * (vy + (k2vy / 2))
    k4vx = step * ax(x + k3x, y + k3y, q1, q2, m)
    k4vy = step * ay(x + k3x, y + k3y, q1, q2, m)
    k4x = step * (vx + k3vx)
    k4y = step * (vy + k3vy)

    # Update the values based on our calculated constants
    vx1[i] = vx = vx + (k1vx + k2vx + k2vx + k3vx + k3vx + k4vx) / 6
    vy1[i] = vy = vy + (k1vx + k2vy + k2vy + k3vy + k3vy + k4vy) / 6
    x1[i] = x = x + ((k1x + 2 * k2x + 2 * k3x + k4x) / 6)
    y1[i] = y = y + ((k1y + 2 * k2y + 2 * k3y + k4y) / 6)

    # Update the time
    time[i] = t = t0 + i * step

Вот функции, которые я использую для ax и ay в предыдущем коде 

def accel_rk4_x(x, y, q1, q2, m):
    const = (q1 * q2) / (4 * math.pi * 8.854e-12 * m)
    return const * (x / ((x ** 2 + y ** 2) ** 1.5))


def accel_rk4_y(x, y, q1, q2, m):
    const = (q1 * q2) / (4 * math.pi * 8.854e-12 * m)
    return const * (y / ((x ** 2 + y ** 2) ** 1.5))

Я ценю любую помощь с этой проблемой! Мне может понадобиться вторая пара глаз.

...