如何在 Julia 中快速求解多项式方程?
How to solve polynomial equations fast in Julia?
我正在研究生物系统的进化模拟。我必须求解多项式方程,找到根 (u*X^3 - N*p*r*X^2 - N*p*X^2 + K^2*u*X - N*K^2*p ), 其中 u 和 K 是常量,N 是常量数组,p, r 是演化参数。基本上,对于每一代人口中的每个人,我需要进行以下计算 (length(N)>>length(p)):
for i = 1:length(p)
for j = 1:length(N)
X[j,i] = mean(fzeros(S -> u*S^3 - p[i]*N[j]*r[i]*S^2 - p[i]*N[j]*S^2 + K^2*u*S - p[i]*K^2*N[j], 0, Inf) )
end
end
我知道我可以通过使用不同内核为不同个体求解方程来并行代码,甚至在每个个体中我也可以并行求解每个 X[j,i]。我想知道处理这种情况的最佳 practice/fast 方法是什么?是否有可能以更快的方式求解单个方程?
这个答案比较了 Roots
、GSL
和 PolynomialRoots
包来解决问题中的问题版本。时间因机器而异,最好在本地 运行。但本质上,GSL
最快,比 Roots
快大约 1000 倍,PolynomialRoots
介于两者之间(但更接近 GSL
)。查看评论了解更多信息。
代码:
# using Julia 0.5.1-pre+4
# Pkg.add("Roots")
# Pkg.add("GSL")
# Pkg.add("PolynomialRoots")
using GSL
using Roots
using PolynomialRoots
function test1(p,r,N,u,K)
X = Matrix{Float64}(length(N),length(p))
for i = 1:length(p)
for j = 1:length(N)
X[j,i] = mean(fzeros(S -> u*S^3 - p[i]*N[j]*r[i]*S^2 - p[i]*N[j]*S^2 + K^2*u*S - p[i]*K^2*N[j], 0, Inf) )
end
end
return X
end
function test2(p,r,N,u,K)
X = Matrix{Float64}(length(N),length(p))
uinv = inv(u)
for i = 1:length(p)
for j = 1:length(N)
X[j,i] = mean(poly_solve_cubic(-uinv*(p[i]*N[j]*r[i]+p[i]*N[j]),K^2,-uinv*p[i]*K^2*N[j]) )
end
end
return X
end
function test3(p,r,N,u,K)
X = Matrix{Float64}(length(N),length(p))
for i = 1:length(p)
for j = 1:length(N)
X[j,i] = mean(filter(x->abs(imag(x))<1.0e-10,
PolynomialRoots.solve_cubic_eq(Complex{Float64}[- p[i]*K^2*N[j], K^2*u, - p[i]*N[j]*r[i] - p[i]*N[j],u])))
end
end
return X
end
K = 1.0
u = 1.0
N = rand(1000)+1.0
p = rand(10)+1.0
r = rand(10)+1.0
res1 = test1(p,r,N,u,K);
res2 = test2(p,r,N,u,K);
res3 = test3(p,r,N,u,K);
using Base.Test
@test_approx_eq res1 res2
@test_approx_eq res1 res3
@time test1(p,r,N,u,K); # Roots
@time test2(p,r,N,u,K); # GSL
@time test3(p,r,N,u,K); # PolynomialRoots
nothing
我的时间:
20.664363 seconds (225.67 M allocations: 13.650 GB, 18.81% gc time) # Roots
0.010303 seconds (80.01 k allocations: 4.044 MB, 75.30% gc time) # GSL
0.215453 seconds (394.90 k allocations: 9.917 MB) # PolynomialRoots
如果您比较这些方法,很高兴听到有关实际问题时间的报告。谢谢。
我正在研究生物系统的进化模拟。我必须求解多项式方程,找到根 (u*X^3 - N*p*r*X^2 - N*p*X^2 + K^2*u*X - N*K^2*p ), 其中 u 和 K 是常量,N 是常量数组,p, r 是演化参数。基本上,对于每一代人口中的每个人,我需要进行以下计算 (length(N)>>length(p)):
for i = 1:length(p)
for j = 1:length(N)
X[j,i] = mean(fzeros(S -> u*S^3 - p[i]*N[j]*r[i]*S^2 - p[i]*N[j]*S^2 + K^2*u*S - p[i]*K^2*N[j], 0, Inf) )
end
end
我知道我可以通过使用不同内核为不同个体求解方程来并行代码,甚至在每个个体中我也可以并行求解每个 X[j,i]。我想知道处理这种情况的最佳 practice/fast 方法是什么?是否有可能以更快的方式求解单个方程?
这个答案比较了 Roots
、GSL
和 PolynomialRoots
包来解决问题中的问题版本。时间因机器而异,最好在本地 运行。但本质上,GSL
最快,比 Roots
快大约 1000 倍,PolynomialRoots
介于两者之间(但更接近 GSL
)。查看评论了解更多信息。
代码:
# using Julia 0.5.1-pre+4
# Pkg.add("Roots")
# Pkg.add("GSL")
# Pkg.add("PolynomialRoots")
using GSL
using Roots
using PolynomialRoots
function test1(p,r,N,u,K)
X = Matrix{Float64}(length(N),length(p))
for i = 1:length(p)
for j = 1:length(N)
X[j,i] = mean(fzeros(S -> u*S^3 - p[i]*N[j]*r[i]*S^2 - p[i]*N[j]*S^2 + K^2*u*S - p[i]*K^2*N[j], 0, Inf) )
end
end
return X
end
function test2(p,r,N,u,K)
X = Matrix{Float64}(length(N),length(p))
uinv = inv(u)
for i = 1:length(p)
for j = 1:length(N)
X[j,i] = mean(poly_solve_cubic(-uinv*(p[i]*N[j]*r[i]+p[i]*N[j]),K^2,-uinv*p[i]*K^2*N[j]) )
end
end
return X
end
function test3(p,r,N,u,K)
X = Matrix{Float64}(length(N),length(p))
for i = 1:length(p)
for j = 1:length(N)
X[j,i] = mean(filter(x->abs(imag(x))<1.0e-10,
PolynomialRoots.solve_cubic_eq(Complex{Float64}[- p[i]*K^2*N[j], K^2*u, - p[i]*N[j]*r[i] - p[i]*N[j],u])))
end
end
return X
end
K = 1.0
u = 1.0
N = rand(1000)+1.0
p = rand(10)+1.0
r = rand(10)+1.0
res1 = test1(p,r,N,u,K);
res2 = test2(p,r,N,u,K);
res3 = test3(p,r,N,u,K);
using Base.Test
@test_approx_eq res1 res2
@test_approx_eq res1 res3
@time test1(p,r,N,u,K); # Roots
@time test2(p,r,N,u,K); # GSL
@time test3(p,r,N,u,K); # PolynomialRoots
nothing
我的时间:
20.664363 seconds (225.67 M allocations: 13.650 GB, 18.81% gc time) # Roots
0.010303 seconds (80.01 k allocations: 4.044 MB, 75.30% gc time) # GSL
0.215453 seconds (394.90 k allocations: 9.917 MB) # PolynomialRoots
如果您比较这些方法,很高兴听到有关实际问题时间的报告。谢谢。