Search code examples
juliasolverode

How to fix "dt <= dtmin. Aborting" error in solveODE

I'm trying to emulate a system of ODEs (Fig3 B in Tilman, 1994.Ecology, Vol.75,No1,pp-2-16) but Julia Integration method failed to give a solution.

The error is dt <= dtmin. Aborting.

using DifferentialEquations
TFour = @ode_def TilmanFour begin
dp1 = c1*p1*(1-p1) - m*p1
dp2 = c2*p2*(1-p1-p2) -m*p2 -c1*p1*p2
dp3 = c3*p3*(1-p1-p2-p3) -m*p3 -c1*p1*p2 -c2*p2*p3
dp4 = c4*p4*(1-p1-p2-p3-p4) -m*p4 -c1*p1*p2 -c2*p2*p3 -c3*p3*p4 
end c1 c2 c3 c4 m

u0 = [0.05,0.05,0.05,0.05]
p = (0.333,3.700,41.150,457.200,0.100)
tspan = (0.0,300.0)
prob = ODEProblem(TFour,u0,tspan,p)
sol = solve(prob,alg_hints=[:stiff])
Solution

  • I think that you read the equations wrong. The last term in the paper is

    sum(c[j]*p[j]*p[i] for j<i)

    Note that every term in the equation for dp[i] has a factor p[i].

    Thus your equations should read 

    dp1 = p1 * (c1*(1-p1) - m)
dp2 = p2 * (c2*(1-p1-p2) - m - c1*p1)
dp3 = p3 * (c3*(1-p1-p2-p3) - m - c1*p1 -c2*p2)
dp4 = p4 * (c4*(1-p1-p2-p3-p4) - m - c1*p1 - c2*p2 - c3*p3)

    where I also made explicit that dpk is a multiple of pk. This is necessary as it ensures that the dynamic stays in the octand of positive variables.

    Using python the plot looks like in the paper

    enter image description here

    def p_ode(p,c,m):
    return [ p[i]*(c[i]*(1-sum(p[j] for j in range(i+1))) - m[i] - sum(c[j]*p[j] for j in range(i))) for i in range(len(p)) ]

c = [0.333,3.700,41.150,457.200]; m=4*[0.100]
u0 = [0.05,0.05,0.05,0.05]
t = np.linspace(0,60,601)

p = odeint(lambda u,t: p_ode(u,c,m), u0, t)
for k in range(4): plt.plot(t,p[:,k], label='$p_%d$'%(k+1));
plt.grid(); plt.legend(); plt.show()