Issue solving ode with solve_ivp
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
Hi, (I mailed this answer twice, the first one seems to have not worked properly...) I believe you made some mistakes in the last 4 integrations. You replace tspan by a range of y (or x) values. What you meant is that you want the solution to stop once it leaves the domain y=[0,4]. However you have an ODE in time, so you cannot know a priori how long it will take for the solution to depart from the saddle point. Here, your y values gets interpreted as a time range (here 4 seconds). When the initial point is very close to the saddle point, the solution departs from it very slowly, so 4 seconds is not sufficient to build the trajectory you are seeking. If you instead ask for a longer integration time (say 5000 s), you get what you want: [image: image.png] To let the simulations stops when the solution leaves the domain of interest, you should use event functions (see the doc of solve_ivp). Here is the last part of you code with my corrections: tspan = np.array([0,5e3]) epsilon = 1e-5 #TODO: ideally, make event functions so that the simulations are stopped when they leave the domain x=[0,4], y=[0,4] to spare computational time ic = [x4-epsilon,y4+epsilon] solution = solve_ivp(f, t_span=tspan, y0=ic, method='RK45', t_eval=None, dense_output=True, args=(a,b,c,d,)) assert solution.success plt.plot(solution.y[0], solution.y[1],'g--', linewidth=4, zorder=80) ic = [x4+epsilon,y4-epsilon] solution = solve_ivp(f, t_span=tspan, y0=ic, method='RK45', t_eval=None, dense_output=True, args=(a,b,c,d,)) assert solution.success plt.plot(solution.y[0], solution.y[1],'m--', linewidth=4, zorder=81) ic = [x4-epsilon,y4-epsilon] solution = solve_ivp(f, t_span=-tspan, y0=ic, method='RK45', t_eval=None, dense_output=True, args=(a,b,c,d,))#, rtol=1e-10) assert solution.success plt.plot(solution.y[0], solution.y[1], 'y', linewidth=4, zorder=82) ic = [x4+epsilon,y4+epsilon] solution = solve_ivp(f, t_span=-tspan, y0=ic, method='RK45', t_eval=None, dense_output=True, args=(a,b,c,d,))#, rtol=1e-10, atol=1E-12) # assert solution.success # this last simulation fails because the point reached is too far from the domain of interest and the solution diverges # however, we already have a nice trajectory inside the required domain plt.plot(solution.y[0], solution.y[1], 'b', linewidth=4, zorder=83) # note that the function g=-f is not needed, simply integrate f in negative time instead Regards, Laurent
Oh my. Thank you so much. This one gave couple of frustrated hours. Also thanks for the other hints. The script was basically a 1:1 translation from the given Matlab file, there however with tspan=[0 100], which didn't work at all, and the MatLab doku did't really explain, how the parameter is used inside their function, it's a black box after all. I think, I understand the solve_ivp better now. Thanks again.
participants (2)
-
Andreas Scholz -
FRANCOIS Laurent