Неожиданные результаты от интегратора scipy.integrate.ode dopri5 (Python) ⇐ Python

[Anonymous]

У меня есть такой скрипт Python:

Код: Выделить всё

def dY_dt(t, Y):
STUFF (defines the dY_dt values)

if dY6_dt < 0.:
print("stuff about current parameters")
1/0 #never fails here
else:
pass

return [dY1_dt, dY2_dt, dY3_dt, dY4_dt, dY5_dt, dY6_dt]

delta_t = 0.01
t_arr = []
Y_arr = []
t_i = 0.00
Y_i = [Y1_i, Y2_i, Y3_i, Y4_i, Y5_i, Y6_i]
t_arr.append(t_i)
Y_arr.append(Y_i)
t_end = 1000.

backend = 'dopri5'
solver = ode(dY_dt).set_integrator(backend)
solver.set_initial_value(Y_i, t_i)

while solver.successful() and (solver.t+delta_t) = 0, for all time, it ought to be the case that dY6dt_total >= 0, for all time.
It is unclear to me why after a full time-step of delta_t = 0.01, that dY6dt_total can be negative. This is for a physical system, and, this happens to correspond to a negative length, which, is non-physical.
The integrator makes a bunch of sub-steps from "t" to "t + delta_t", and for one crash of the code, the crash occurs after the time-step from t = 84.01 to t = 84.02. All of the values of dY6_dt during each of these sub-steps are positive, as they should be, but, that suggests the integrator isn't behaving as I would expect (rather than an error in defining dY6_dt, causing it to be negative).
Here are those values:
[code]t = 84.01
Y6 = 20160.4702157939
dY6_daddet = 2.499582626169507e-05

t = 84.012
Y6 = 20160.470215843894
dY6_dt = 3.1744418995557098

t = 84.013
Y6 = 20160.47735830693
dY6_dt = 2.8810884888089463e-05

t = 84.018
Y6 = 20160.35170456511
dY6_dt = 0.0011118447521629496

t = 84.0188889
Y6 = 20160.10211440677
dY6_dt = 3.2265511820351254

t = 84.02
Y6 = 20160.11990354711
dY6_dt = 3.2239658437410013

t = 84.02
Y6 = 20160.4640434123
dY6_dt = 3.1749812457108955

But, after the integration for this time-step:

Код: Выделить всё

dY6dt_total = -0.006172381603391841

Running this script with the same parameters and initial conditions, but, with delta_t = 0.001, and, delta_t = 0.0001, did not prevent "dY6dt_total < 0." from occurring around t ~ 80-90, it just changed the exact time it occurred.
Can anyone help with this? Is there something I've misunderstood about how 'dopri5' works? I cannot see how if dY6_dt >= 0. for all time, that after a time-step the value of Y6 decreases.


Источник: https://stackoverflow.com/questions/650 ... ri5-python