Hi everyone,
I have an issue with integrate.solve_ivp. In a Matlab example with the Matlabs function
[~,x] = ode45(@(t,x) g(t,x,a,b,c,d),tspan,ic), I am able to create this output, visible the 4 solution lines from the saddle point.
[cid:87161cf7-f5c2-47d6-a3a6-75c4b17157ce]
The problematic parts are the last 4 blocks calculating solutions from the saddle point at [1,1], the unstable node [0,0] and the stable nodes [0,c] and [a,0]. This function is defined as follows:
def lotkavolterra_variant(a=3, b=2, c=2, d=1, xmin=0, xmax=4, ymin=0, ymax=4, stepsize=.25):
## Lotka-Volterra system with intraspecific competition
# Skript Figure 4.27
x = np.arange(xmin, xmax+stepsize, stepsize)
y = np.arange(ymin, ymax+stepsize, stepsize)
yval = np.arange(ymin, ymax+stepsize*0.1, stepsize*0.1) # evaluate curve 4 times more often to get smoothness
[X, Y] = np.meshgrid(x, y)
dx = X*(a - X - b*Y)
dy = Y*(c - Y - d*X)
f = lambda t, x,a,b,c,d : [x[0]*(a - x[0] - b*x[1]), x[1]*(c - x[1] - d*x[0])]
g = lambda t, x,a,b,c,d : [-x[0]*(a - x[0] - b*x[1]), -x[1]*(c - x[1] - d*x[0])]
plt.figure(figsize = (13,13))
plt.quiver(X, Y, dx, dy, angles = 'xy', color='#0086b3', width=0.0015)
plt.grid()
plt.xlabel('x')
plt.ylabel('y')
plt.axis([xmin,xmax,ymin,ymax])
tspan=[0,50]
x0stepsize = .0001
for x0 in np.arange(0, .01+x0stepsize, x0stepsize):
ic = [x0, .01-x0]
solution = solve_ivp(f, [xmin, xmax], ic, method='RK45', t_eval=yval, dense_output=True, args=(a,b,c,d,))
plt.plot(solution.y[0], solution.y[1],'r')
x0stepsize = .1
for x0 in np.arange(0, 10+x0stepsize, x0stepsize):
ic = [x0, 10-x0]
solution = solve_ivp(f, [xmin, xmax], ic, method='RK45', t_eval=yval, dense_output=True, args=(a,b,c,d,))
plt.plot(solution.y[0], solution.y[1],'r')
x4 = (a - b*c)/(1 - b*d)
y4 = (c - a*d)/(1 - b*d)
plt.plot( 0, 0 ,'go',fillstyle='none', ms=15, zorder=100)[0].set_clip_on(False)
plt.plot (0, c ,'yo',fillstyle='full', ms=15, zorder=100)[0].set_clip_on(False)
plt.plot( a, 0 ,'mo',fillstyle='full', ms=15, zorder=100)[0].set_clip_on(False)
plt.plot(x4,y4 ,'bo',fillstyle='none', ms=15, zorder=100)[0].set_clip_on(False)
epsilon = 0.00001
ic = [x4-epsilon,y4+epsilon]
solution = solve_ivp(f, [xmin, xmax], ic, method='RK45', t_eval=yval, dense_output=True, args=(a,b,c,d,))
plt.plot(solution.y[0], solution.y[1],'g--', linewidth=2, zorder=80)
ic = [x4+epsilon,y4-epsilon]
solution = solve_ivp(f, [xmin, xmax], ic, method='RK45', t_eval=yval, dense_output=True, args=(a,b,c,d,))
plt.plot(solution.y[0], solution.y[1],'m--', linewidth=2, zorder=81)
ic = [x4-epsilon,y4-epsilon]
solution = solve_ivp(g, [ymin, ymax], ic, method='RK45', t_eval=yval, dense_output=True, args=(a,b,c,d,))
plt.plot(solution.y[0], solution.y[1], 'y', linewidth=2, zorder=82)
ic = [x4+epsilon,y4+epsilon]
solution = solve_ivp(g, [xmin, xmax], ic, method='RK45', t_eval=yval, dense_output=True, args=(a,b,c,d,))
plt.plot(solution.y[0], solution.y[1], 'b', linewidth=2, zorder=83)
While Matlab is able to compute them, the Scipy implementation failed. I tried lots of variations of methods, first_step, max_step (I even tried to set max_step to 0.0001), rtol and atol adjustments. The result is always like the image below. It seems the computation stops after about 40 steps and a short distance
[cid:e6b716a4-b876-4136-9b84-281d3810d8c8]
unless I set epsilon to something quite high like 0.1. Maybe I'm missing something, but I think something is not right with the Scipy implementation. The same issue also arises with other examples like a pendulum calculation which should calculate the whole definition space but stops midair.
[cid:bec1cd4a-d8b4-422f-acfc-d24ed446a147]
Thanks in advance and kind regards