求解预定义函数的非线性常微分方程西皮.奥丁

Question

我想解一个非线性常微分方程

θ2=（C+j（θ2））**-1*（f（t）–g（θ1）–h（θ0））

其中f（）、g（）、h（）和j（）是已经定义的函数，它们将θ2、θ1、θ0或t作为输入。θ2和θ1是θ0随时间t的二阶和一阶导数

我一直在用SciPy.odeint公司函数使用以下代码：

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)

Answer 1

方程式可改写为：

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

其中a和b是常数，F对应于（F（t）–g（θ1）–h（θ0））。你知道吗

它是θ“”的二阶多项式函数，有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();