Использование Рунге-Кутты 4-го порядка для решения дифференциального уравнения 2-го порядка затухающего осциллятора

Я думаю, что может возникнуть некоторая путаница с термином внешнего воздействия и производной вспомогательной функцией Рунге-Кутты F. F в RK4 возвращает производную dX/dt системы дифференциальных уравнений первого порядка X. Вынуждающий член в демпфированном осцилляторе, к сожалению, также называется F, но это функция t.

Одна из ваших проблем заключается в том, что арность (количество параметров) вашей функции RungeKutta() и ваш вызов этой функции не совпадают: вы пытались сделать RungeKutta(F, y0, t, func), но функция RungeKutta() принимает аргументы (f, y0, x) только в таком порядке.

Другими словами, параметр f в вашей текущей RungeKutta() функции должен инкапсулировать принудительную функцию F(t).

Вы можете сделать это с помощниками:

# A constant function in your case, but this can be any function of `t`
def applied_force(t):
    # Note, you did not provide a value for `omega0`
    return omega0 ** 2

def rk_derivative_factory(c, k, F):
    return lambda X, t: np.array([X[1], -c * X[1] - k * X[0] + F(t)])

rk_derivative_factory() — это функция, которая принимает коэффициент демпфирования, жесткость пружины и функцию принуждения F(t), а возвращает функцию принимает в качестве аргументов систему X и временной шаг t (поскольку это то, что требуется от вас реализацией RungeKutta()) .

Затем,

omega0 = 0.234
m, c, k = 1, 0.25, 2
oscillator = HarmonicOscillator(m, c, k)

f = rk_derivative_factory(oscillator, applied_force)
x_osc = oscillator.solve(f, [1, 0])

Где solve() определяется так:

def solve(self, func, y0):
    t = np.linspace(0,10 * self.period(), 1000)
    return RungeKutta(f, y0, t)

Кроме того, я настоятельно рекомендую быть более последовательным в именах переменных. Вы назвали начальное состояние вашего осциллятора x0 и передавали его RungeKutta() в качестве аргумента для параметра y0, а параметр xRungeKutta() представляет время... Довольно запутанно.

Полное решение

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

Обратите внимание, что одна вещь, которую вы, возможно, захотите рассмотреть, — заставить функцию HarmonicOscillator.solve() принимать решатель. Затем вы можете поиграть с разными интеграторами.

import numpy as np


def RungeKutta(f, y0, x):
    y = np.zeros((len(x), len(y0)))
    y[0, :] = np.array(y0)
    h = x[1] - x[0]
    for i in range(0, len(x) - 1):
        # Many slight changes below
        k1 = np.array(f(y[i, :], x[i]))
        k2 = np.array(f(y[i, :] + h * k1 / 2, x[i] + h / 2))
        k3 = np.array(f(y[i, :] + h * k2 / 2, x[i] + h / 2))
        k4 = np.array(f(y[i, :] + h * k3, x[i] + h))
        y[i + 1, :] = y[i, :] + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4)
    return y


# A constant function in your case, but this can be any function of `t`
def applied_force(t):
    # Note, you did not provide a value for `omega0`
    return omega0 ** 2


def rk_derivative_factory(osc, F):
    return lambda X, t: np.array([X[1], (F(t) - osc.c * X[1] - osc.k * X[0]) / osc.m])


class HarmonicOscillator:
    def __init__(self, m, c, k):
        if (type(m) in (int, float)) and (m > 0):
            self.m = m 
        else:
            raise ValueError("Parameter 'm' must be a positive number")
        if (type(c) in (int, float)) and (c > 0):
            self.c = c
        else:
            raise ValueError("Parameter 'c' must be a positive number")
        if (type(k) in (int, float)) and (k > 0):
            self.k = k 
        else:
            raise ValueError("Parameter 'k' must be a positive number")
        
        self.T = 2 * np.pi * (self.m / self.k)**(0.5)

    def period(self):
        return self.T

    def solve(self, func, y0):
        t = np.linspace(0, 10 * self.period(), 1000)
        return RungeKutta(func, y0, t)

Демо:

import plotly.graph_objects as go

omega0 = 0.234
m, c, k = 1, 0.25, 2
oscillator = HarmonicOscillator(m, c, k)

f = rk_derivative_factory(oscillator, applied_force)
x_osc = oscillator.solve(f, [1, 0])

x, dx = x_osc.T
t = np.linspace(0, 10 * oscillator.period(), 1000)

fig = go.Figure(go.Scatter(x=t, y=x, name = "x(t)"))
fig.add_trace(go.Scatter(x=t, y=dx, name = "x'(t)"))

Выход: