使用 scipy.odeint 求解预定义函数的非线性 ODE答案

【问题标题】：Solve Non linear ODE of predefined functions with scipy.odeint使用 scipy.odeint 求解预定义函数的非线性 ODE
【发布时间】：2020-01-14 02:26:45
【问题描述】：

我想求解形式为

的非线性常微分方程

Theta2 = (C + j(Theta2))**-1 * (f(t) – g(Theta1) -h(Theta0))

其中 f()、g()、h() 和 j() 是已经定义的函数，它们以 Theta2、Theta1、Theta0 或 t 作为输入。 Theta2 和 Theta1 是 Theta0 随时间 t 的二阶和一阶导数。

我一直在使用以下代码使用 SciPy.odeint 函数求解没有 j(Theta2) 项的方程：

from scipy.integrate import odeint

def ODE():
    def g(Theta, t):
        Theta0 = Theta[0]
        Theta1 = Theta[1]
        Theta2 = (1/C)*( f(t) - g(Theta1)  - h(Theta0))

        return Theta1, Theta2
    init = 0, 0 # Initial conditions on theta0 and theta1 (velocity) at t=0
    sol=odeint(g, init, t)
    A = sol[:,1]
    B = sol[:,0]
    return(A, B)

【问题讨论】：

这里是一个起点：stackoverflow.com/a/50796835/8069403你能提供更多关于函数j(Theta2)的信息吗？
@xdze2 J(Theta2) 仅将 Theta2 作为输入，并为 Theta 2 的每个值生成一个输出。该函数由一系列常数乘以 Theta2 组成。简化了它的 C * Theta2 = 输出。我已经编辑了缺少 ^-1 的代码，这会使问题变得不清楚。

标签： python scipy ode

【解决方案1】：

方程可以改写为：

          F(t, theta, theta')
theta'' = -------------------
            a + b*theta''

其中 a 和 b 是常数，F 对应于 (f(t) – g(Theta1) -h(Theta0))。

它是 theta'' 的二阶多项式函数，有 2 个解（考虑 b!=0 和 a^2 + 4*b*F>0）：

theta'' = -( sqrt(a^2 + 4*b*F) +/- a )/(2*b)

这个新方程的形式为 y' = f(t, y)，可以使用常规 ODE 求解器求解。

这是一个使用solve_ivp 的示例，它是odeint 的替代品：

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

a = 20
b = 1

def f(t, y, dydt):
    return t + y**2

def ode_function_plus(t, Y):
    y = Y[0]
    dydt = Y[1]

    d2y_dt2 = -(np.sqrt(a**2 + 4*b*f(t, y, dydt)) + a )/(2*b)
    return [dydt, d2y_dt2]

def ode_function_minus(t, Y):
    y = Y[0]
    dydt = Y[1]

    d2y_dt2 = -(np.sqrt(a**2 + 4*b*f(t, y, dydt)) - a )/(2*b)
    return [dydt, d2y_dt2]

# Solve
t_span = [0, 4]
Y0 = [10, 1]
sol_plus = solve_ivp(ode_function_plus, t_span, Y0)
sol_minus = solve_ivp(ode_function_minus, t_span, Y0)
print(sol_plus.message)

# Graph
plt.plot(sol_plus.t, sol_plus.y[0, :], label='solution +a');
plt.plot(sol_minus.t, sol_minus.y[0, :], label='solution -a');
plt.xlabel('time'); plt.ylabel('y'); plt.legend();

【讨论】：