Я решаю задачу с помощью RK4 в Python. Я почти решил, но результаты не совпадают. Может ли кто-нибудь, пожалуйста, мне помочь?

это DOI статьи: 10.1017/S0022112003003835, система находится на странице 190, а процедура на странице 191.

import numpy as np
import matplotlib.pyplot as plt

# Параметры
R = 1000                  # Число Рейнольдса
Delta = 0.3               # Амплитуда колебаний
L = 1.0                   # Длина области
N = 101                  # Количество узлов сетки
eta = np.linspace(0, L, N)  # Пространственная сетка
h = eta[1] - eta[0]       # Шаг сетки
T_end = 10000.0            # Общее время моделирования
num_time_steps = 10000      # Количество шагов по времени (увеличено для более точного разрешения)
dt = T_end / num_time_steps  # Шаг по времени

# Функции движения стенки
H = lambda t: 1 + Delta * np.cos(2 * t)
H_dot = lambda t: -2 * Delta * np.sin(2 * t)

# Инициализация F и G
F = np.zeros(N)
G = np.zeros(N)

# Вспомогательные функции для пространственных производных
def dFdx(F):
    """Вычислить первую производную с использованием центральных разностей."""
    dF = np.zeros_like(F)
    dF[1:-1] = (F[2:] - F[:-2]) / (2 * h)
    return dF

def d2Fdx2(F):
    """Вычислить вторую производную с использованием центральных разностей."""
    d2F = np.zeros_like(F)
    d2F[1:-1] = (F[2:] - 2 * F[1:-1] + F[:-2]) / h**2
    return d2F

def spatial_operator_G(G):
    """Применить оператор L к G."""
    d2G = d2Fdx2(G)
    dG = dFdx(G)
    L_G = d2G + dG / eta - G / eta**2
    return L_G

# Правая часть для Gt
def compute_rhs(F, G, t):
    """Вычислить правую часть для F и G."""
    dF = dFdx(F)
    d2F = d2Fdx2(F)
    dG = dFdx(G)
    L_G = spatial_operator_G(G)

    RHS_F = d2F + dF / eta - F / eta**2 + H(t)**2 * G
    RHS_G = (
        -H_dot(t) / H(t) * eta * dG
        - F * dG / H(t)
        + dF * G / H(t)
        + 2 * F * G / (H(t) * eta)
        + L_G / (H(t)**2 * R)
    )
    return RHS_F, RHS_G

# Интеграция методом Рунге-Кутты 4-го порядка
F_sol = []  # Хранение решений F
G_sol = []  # Хранение решений G
times = []  # Хранение шагов времени

t = 0.0
for n in range(num_time_steps):
    F_sol.append(F.copy())
    G_sol.append(G.copy())
    times.append(t)

    # Шаги метода Рунге-Кутты
    k1_F, k1_G = compute_rhs(F, G, t)
    k2_F, k2_G = compute_rhs(F + 0.5 * dt * k1_F, G + 0.5 * dt * k1_G, t + 0.5 * dt)
    k3_F, k3_G = compute_rhs(F + 0.5 * dt * k2_F, G + 0.5 * dt * k2_G, t + 0.5 * dt)
    k4_F, k4_G = compute_rhs(F + dt * k3_F, G + dt * k3_G, t + dt)

    F += dt * (k1_F + 2 * k2_F + 2 * k3_F + k4_F) / 6
    G += dt * (k1_G + 2 * k2_G + 2 * k3_G + k4_G) / 6

    # Граничные условия
    F[0], G[0] = 0, 0                  # Условие Дирихле при eta = 0
    F[-1] = -H_dot(t + dt)             # Условие Дирихле при eta = 1
    F[-2] = F[-1] - h * H(t + dt)      # Условие Неймана при eta = 1

    # Инкремент времени
    t += dt

# Преобразование решений в массивы
F_sol = np.array(F_sol)
G_sol = np.array(G_sol)
times = np.array(times)

# Построение результатов для временного диапазона от t = 8100 до t = 8200
time_range = np.where((times >= 8100) & (times <= 8200))

plt.figure(figsize=(10, 6))
plt.plot(times[time_range], np.squeeze(F_sol[time_range, -1]), label="F при eta=1/4")
plt.plot(times[time_range], np.squeeze(G_sol[time_range, -1]), label="G при eta=1/4")
plt.xlabel('Время t')
plt.ylabel('Значения решений')
plt.title('Эволюция F и G во времени от t = 8100 до t = 8200 при eta=1/4')
plt.legend()
plt.grid()
plt.show()