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
