Как использовать сгенерированную Якоби матрицу в решении системы ОДУ? - PullRequest
Новая
       13

Как использовать сгенерированную Якоби матрицу в решении системы ОДУ?

0 голосов
/ 02 января 2019

У меня есть система ODE первого порядка, которая состоит из 3 разн.уравнение х.Я хочу решить эту проблему с помощью метода BDF scipy.integrate.solve_ivp.Поэтому мне нужно рассчитать матрицу Якоби системы (и сделать это с помощью SymPy).

Если я не понял неправильно;согласно документу scipy.integrate.solve_ivp, вы должны ввести матрицу jacobien в виде jac (t, u), где u должно быть переменными состояния вашей системы ODE.С этой целью я правильно раскладываю матрицу Якобиана.

И здесь возникает моя проблема.Хотя я могу вычислить jac (t, u) с некоторыми (t, u), такими как ((1/800), (150,1E-6,3)), я не смог отправить аргументы массива в мой jac.когда я представляю jac (t, u) в качестве аргумента для solve_ivp, выдается сообщение об ошибкеТак как же мне ввести матрицу JAC?Или мой lambdify не подходит?

Это мой код.Любая помощь, я ценю это.

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

def cvs(t,u):
    u1,u2,u3 = u
    def Qmi(t):
        return t**2
    u1p = Qmi(t)*u3
    u2p = (u1**2)*np.cos(2*np.pi*200*t)
    u3p = (np.sin(2*np.pi*t))*u2**-1
    return [u1p,u2p,u3p]
def jac_func():
    ######### DEFINE THE ODE SYSTEM #########
    import sympy
    sympy.init_printing()
    t = sympy.symbols("t")
    Q_mi = sympy.Function("Q_mi")(t)
    u1 = sympy.Function("u1")(t)
    u2 = sympy.Function("u2")(t)
    u3 = sympy.Function("u3")(t)
    Q_mi = t**2
    u1p = (u3*Q_mi)
    u2p = (u1**2)*sympy.cos(2*sympy.pi*200*t)
    u3p = sympy.sin(2*sympy.pi*5*t)*u2**-1
    ####### CALCULATE THE JACOBIEN ########
    ode_rhs = sympy.Matrix([u1p,u2p,u3p])
    ode_var = sympy.Matrix([u1,u2,u3])
    jac = sympy.Matrix([[ode.diff(var) for var in ode_var]for ode in ode_rhs])
    u = (u1,u2,u3)
    jac_np = sympy.lambdify((t,u),jac,"numpy")
    return jac_np

jac_np = jac_func()
U_0 = [500,20,20]
t = np.linspace(0,100,10000)

solf = solve_ivp(cvs,(0,100),y0=U_0,method = 'BDF',jac=jac_np(t,U_0),t_eval=t)

сообщение об ошибке: 

ValueError                                Traceback (most recent call last)
<ipython-input-1-8b86ffb3a7cf> in <module>()
41 t = np.linspace(0,100,10000)
42 
---> 43 solf = solve_ivp(cvs,(0,100),y0=U_0,method = 'BDF',jac=jac_np(t,U_0),t_eval=t)

<lambdifygenerated-1> in _lambdifygenerated(t, _Dummy_188)
  1 def _lambdifygenerated(t, _Dummy_188):
  2     [_Dummy_185, _Dummy_186, _Dummy_187] = _Dummy_188
----> 3     return (array([[0, 0, t**2], [2*_Dummy_185*cos(400*pi*t), 0, 0], [0, -sin(10*pi*t)/_Dummy_186**2, 0]]))

ValueError: setting an array element with a sequence.

1 Ответ

0 голосов
/ 04 января 2019

Вы получаете проблему, потому что вы делаете то, что говорится в сообщении об ошибке, вы передаете массив, где процедура ожидает одно число. В 

solf = solve_ivp(cvs,(0,100),y0=U_0,method = 'BDF',jac=jac_np(t,U_0),t_eval=t)

вы пытаетесь использовать постоянную матрицу jac_np(t,U_0) для аргумента якобиана. Однако в этот момент t содержит все значения t, из которых вы хотите получить сэмплы. Список [ array, scalar, scalar ] не совместим с массивами numpy.

Короче говоря, уберите аргументы, передайте якобиан как вызываемую функцию, как вы, вероятно, и предполагали, 

solf = solve_ivp(cvs,(0,100),y0=U_0,method = 'BDF',jac=jac_np, t_eval=t)
...