Значение погрешности измерения при использовании решения Scipy Ode - PullRequest
0 голосов
/ 19 апреля 2020

В настоящее время пытаюсь запустить следующий код

import numpy as np

from scipy.interpolate import interp1d

from scipy.integrate import ode


def dydt(t,y,tsi,rho):

    Lambda      = 10^-4
    beta        = 0.0065
    lambda2      = 0.08

    rho = interp1d(tsi,rho, fill_value = 'extrapolate')


    #one group delayed precursor
    dydt1 = (-lambda2*y[0] + beta*y[1])
    #power
    dydt2 = (((rho(t)-beta)/Lambda)*y[1]+(lambda2*y[0])/Lambda)
    return dydt1, dydt2


x = 21

dt = 1

tsi=np.arange(0,x,dt)
dt = [1]

rho=np.ones(x)*0.0025
y0= [1,0]
t0 = [0,x-1]


sol = []

i = ode(dydt)

i.set_integrator('dopri5')

i.set_initial_value(t0,y0)

i.set_f_params(tsi,rho)


for t in tsi[1:]:

    i.integrate(i.t+dt)

    sol.append(i.y)

Но я получаю эту ошибку:

runfile ('C: / Users / --- - / Desktop / ---. Py ', wdir =' C: / Users / --- / Desktop ') Traceback (последний последний вызов):

File "C: \ Users --- \ Desktop---.py ", строка 46, в i.integrate (i.t + dt)

Файл" C: \ Users --- \ anaconda3 \ lib \ site- пакеты \ scipy \ integrate_ode.py ", строка 432, в файле интегрировать self.f_params, self.jac_params)

Файл" C: \ Users --- \ anaconda3 \ lib \ site-packages \ scipy \ integrate_ode.py ", строка 1172, в кортеже выполнения (self.call_args) + (f_params,)))

ValueError: 0-е измерение должно быть 2, но получено 0 (не определено).

Я следовал аналогичному примеру онлайн, но главное отличие в том, что я возвращаю 2 уравнения в моей функции, и мои два аргумента являются массивами. Я пытался использовать odeint и solve_ivp, но они дают крайне неточные ответы.

Спасибо.

1 Ответ

0 голосов
/ 19 апреля 2020

Еще раз проверьте строку

i.set_initial_value(t0,y0)

в документации, аргумент времени идет вторым (потому что это необязательный аргумент со значением по умолчанию 0). И начальное время - скаляр, а не кортеж.

Вы должны определить функцию интерполяции только один раз, а затем передать только это как параметр или даже оставить его как глобальный объект.

sol.append(i.y)

может также не выполнять то, что вы думаете, в векторном значении он создает массив идентичных копий ссылок на i.y. Используйте i.y.copy(), чтобы получить копию, которая удалена из внутренних частей решателя.

...