Использование решите_ivp вместо odeint, чтобы решить начальное значение проблемы - PullRequest
Новая
       79

Использование решите_ivp вместо odeint, чтобы решить начальное значение проблемы

0 голосов
/ 21 ноября 2018

В настоящее время я решаю следующую систему уравнений ODE, используя odeint

dx / dt = (-x + u) /2.0

dy / dt = (-y + x) /5.0

начальные условия: x = 0, y = 0

Тем не менее, я хотел бы использовать solve_ivp, который кажется рекомендуемым вариантом для этого типа проблем, но, честно говоря, я незнаете, как адаптировать код ...

Вот код, который я использую с odeint: 

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

def model(z, t, u):
    x = z[0]
    y = z[1]
    dxdt = (-x + u)/2.0
    dydt = (-y + x)/5.0
    dzdt = [dxdt, dydt]
    return dzdt

def main():
    # initial condition
    z0 = [0, 0]

    # number of time points
    n = 401

    # time points
    t = np.linspace(0, 40, n)

    # step input
    u = np.zeros(n)
    # change to 2.0 at time = 5.0
    u[51:] = 2.0

    # store solution
    x = np.empty_like(t)
    y = np.empty_like(t)
    # record initial conditions
    x[0] = z0[0]
    y[0] = z0[1]

    # solve ODE
    for i in range(1, n):
        # span for next time step
        tspan = [t[i-1], t[i]]
        # solve for next step
        z = odeint(model, z0, tspan, args=(u[i],))
        # store solution for plotting
        x[i] = z[1][0]
        y[i] = z[1][1]
        # next initial condition
        z0 = z[1]

    # plot results
    plt.plot(t,u,'g:',label='u(t)')
    plt.plot(t,x,'b-',label='x(t)')
    plt.plot(t,y,'r--',label='y(t)')
    plt.ylabel('values')
    plt.xlabel('time')
    plt.legend(loc='best')
    plt.show()

main()

1 Ответ

0 голосов
/ 21 ноября 2018

Важно, чтобы solve_ivp ожидал f (t, z) в качестве правой части ODE.Если вы не хотите менять свою функцию ode, а также хотите передать свой параметр u, я рекомендую определить функцию-обертку: 

def model(z, t, u):
    x = z[0]
    y = z[1]
    dxdt = (-x + u)/2.0
    dydt = (-y + x)/5.0
    dzdt = [dxdt, dydt]
    return dzdt

def odefun(t, z):
    if t < 5:
        return model(z, t, 0)
    else:
        return model(z, t, 2)

Теперь легко вызывать solve_ivp: 

def main():
    # initial condition
    z0 = [0, 0]

    # number of time points
    n = 401

    # time points
    t = np.linspace(0, 40, n)

    # step input
    u = np.zeros(n)
    # change to 2.0 at time = 5.0
    u[51:] = 2.0

    res = solve_ivp(fun=odefun, t_span=[0, 40], y0=z0, t_eval=t)
    x = res.y[0, :]
    y = res.y[1, :]

    # plot results
    plt.plot(t,u,'g:',label='u(t)')
    plt.plot(t,x,'b-',label='x(t)')
    plt.plot(t,y,'r--',label='y(t)')
    plt.ylabel('values')
    plt.xlabel('time')
    plt.legend(loc='best')
    plt.show()

main()

Обратите внимание, что без передачи t_eval=t решатель автоматически выберет моменты времени внутри tspan, в которые будет сохранено решение.

...