Как исправить ошибку «dt <= dtmin. Aborting» в solveODE - PullRequest
Новая
       65

Как исправить ошибку «dt <= dtmin. Aborting» в solveODE

0 голосов
/ 01 февраля 2019

Я пытаюсь эмулировать систему ODE (Fig3 B в Tilman, 1994. Экология, Vol.75, No1, pp-2-16), но метод интеграции Джулии не дал решения.

Ошибка dt <= dtmin.Отмена. </p> 

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])

1 Ответ

0 голосов
/ 02 февраля 2019

Я думаю, что вы неправильно прочитали уравнения.Последний член в статье равен 

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

. Обратите внимание, что каждый член в уравнении для dp[i] имеет коэффициент p[i].

Таким образом, ваши уравнения должныпрочитайте 

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)

, где я также указал, что dpk - это кратное pk.Это необходимо, поскольку гарантирует, что динамика остается в октане положительных переменных.

При использовании python график выглядит так, как в 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()
...