0

I would like to use solve_ivp for simulating a system which evolves the following differential function:

dy(t) = AD(t) + By(t)

As for me, t is integer counted from 0, I did

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

def fun(t, y, D, A, B):
    dD1 = A*D[int(t)] - B*y
    return dD1#[dD1, dD2]

D = np.zeros((1024))
D[10:34] = 12
tspan = np.linspace(0, np.size(D)-1, np.size(D))

sol = solve_ivp(lambda t, y: fun(t, y, D, A=0.8, B=0.025),
                [tspan[0], tspan[-1]], 
                [0], method='RK45', t_eval=tspan)

plt.figure()
plt.plot(sol.t, sol.y[0,])

It gives me the result expected with a peak of y at t=34. enter image description here However, if I shift the interval at which D has non-zero values by simply replacing

D[10:34] = 12

with

D[420:444] = 12

I expect that I can get a y of the same shape but only shifted, while the result gives me all 0... In addition, sometimes I meet the warning message:

RuntimeWarning: divide by zero encountered in double_scalars max(1, SAFETY * error_norm ** (-1 / (order + 1))))

Fay
  • 105
  • 1
  • 5

1 Answers1

0

I guess that I've found the reason, simply caused by step size, so pass the parameters

first_step=np.size(D)/100, max_step=np.size(D)/100

to solve_ivp function solves the problem.

Fay
  • 105
  • 1
  • 5