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.

 



 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 

 

 


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.



 

Thanks in advance and kind regards