时间戳

首页 > 解决方案 > 使用 scipy.odeint 求解预定义函数的非线性 ODE

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

问题描述

我想求解形式为的非线性常微分方程

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)

标签： pythonscipyode

解决方案

方程可以重写为：

          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();

推荐阅读