使用 Runge-Kutta 方法和牛顿方法求解 ODEs

Question

大家。我正在尝试求解一个初始条件未知的方程式。基本上我想做的是求解图中给出的常微分方程。

提出这个方程的作者提供的计算方法：用4阶Runge-Kutta方法求解。根据f′(∞)=0,θ(∞)=0,φ(∞)=0(constraints)，未知的初始条件可以通过牛顿法得到。

在我的代码中我把图中的常量写成了对应的数字。代码中的

system_of_odes

就是问题最终需要求解的方程式。对应的初始条件为x(1~8) = [0, 0, u[0], -4, 1, u[1], 1, u[2]], 其中u[0] u[1] u [2] 是需要根据约束确定的未知数。

我尝试应用牛顿法找到最符合约束条件（f′（∞）=0，θ（∞）=0，φ（∞）=0）的初始条件（求解 u）。

我写的代码如下：

import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import differential_evolution
import time
from extensisq import BS5, SWAG
from multiprocessing import Process, Manager

def system_of_odes(t, x, u):
    x1, x2, x3, x4, x5, x6, x7, x8 = x
    u0, u1, u2 = u
    dx1 = 1
    dx2 = x3
    dx3 = x4
    dx4 = -x4 * x2 + 10 * x3 + (x3) ** 2
    dx5 = x6
    dx6 = -10 * (x2 * x6 - 2 * x3 * x5) - 0.001 * (x6) ** 2 - x6 * x8
    dx7 = x8
    dx8 = -1 * (x2 * x8 - 2 * x3 * x7) - 0.001 / 1 * ( -10 * (x2 * x6 - 2 * x3 * x5) - 0.001 * (x6) ** 2 - x6 * x8)
    return [dx1, dx2, dx3, dx4, dx5, dx6, dx7, dx8]

def solve_ivp_with_timeout(u, result_dict):
    x0 = [0.0, 0.0, u[0], -2.0 * (1.0 + 1.0), 1.0, u[1], 1.0, u[2]]
    sol = solve_ivp(system_of_odes, (0, 10), x0, args=(u,), method='LSODA')
    result_dict['sol'] = sol

def objective(u, result_dict=None):
    manager = Manager()
    result_dict = manager.dict()
    p = Process(target=solve_ivp_with_timeout, args=(u, result_dict))
    p.start()
    p.join(timeout=30)

    if p.is_alive():
        p.terminate()
        penalty = 1e6
        print('terminate')
        return np.array([1e6, 1e6, 1e6])

    sol = result_dict['sol']
    x3_inf = sol.y[2, -1]
    x5_inf = sol.y[4, -1]
    x7_inf = sol.y[6, -1]

    return np.array([x3_inf, x5_inf, x7_inf])

def jacobian(u):
    eps = 1e-5
    J = np.zeros((3, 3))
    for i in range(3):
        u_plus = u.copy()
        u_plus[i] += eps
        u_minus = u.copy()
        u_minus[i] -= eps
        J[:, i] = (objective(u_plus) - objective(u_minus)) / (2 * eps)
    return J

def newton_method(u0, max_iter=300, tol=1e-3):
    u = u0.copy()
    delta_u = None

    print(f"Initial guess: u = {u0}")

    for i in range(max_iter):
        f = objective(u)
        J = jacobian(u)
        delta_u = np.linalg.pinv(J) @ (-f)
        u += delta_u
        print(f"Iteration {i + 1}: u = {u}, delta_u = {delta_u}, objective = {np.linalg.norm(f)}")
        if np.linalg.norm(delta_u) < tol:
            break

    if np.linalg.norm(delta_u) >= tol:
        print("Warning: not converge.")
    return u

def fitness(u):
    start_time = time.time()
    u0 = newton_method(u)
    end_time = time.time()
    runtime = end_time - start_time
    obj_value = np.linalg.norm(objective(u))

    penalty = 0
    if obj_value > 1:
        penalty += 1000 * obj_value
    if runtime > 30:
        penalty += 100 * (runtime - 1)

    return obj_value + penalty

def main():
    bounds = [(0, 2), (-200, 200), (-200, 200)]
    result = differential_evolution(fitness, bounds)
    u0 = result.x
    u = newton_method(u0)
    print("Solution:")
    print(f"u = {u}")

if __name__ == "__main__":
    main()

当我在solve ivp中使用LSODA时，会出现以下警告：

lsoda--  warning..internal t (=r1) and h (=r2) are
       such that in the machine, t + h = t on the next step  
       (h = step size). solver will continue anyway
      in above,  r1 =  0.1052021784050D+01   r2 =  0.1007877569584D-15

和

 UserWarning: lsoda: Repeated convergence failures (perhaps bad Jacobian or tolerances).
  warnings.warn('{:s}: {:s}'.format(self.__class__.__name__,

同时，除了尝试使用losda方法外，我还尝试使用BS5、SWAG、Radau、RK45等各种求解ivp的方法

SWAG方法的结果如下图所示，所有的方法都不能正确解决这个问题

我也试过 MATLAB，但效果也不太好。

希望大家能指出我的错误或提出一些建议，帮助我在约束条件下正确求解u值和ODEs。

非常感谢您的帮助。