julia - 如何在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是不断发展的参数。基本上，对于每一代人口中的每个人，我都需要进行以下计算（长度（N）>>长度（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]。我想知道解决这种情况的最佳实践/快速方法是什么？是否可以以更快的方式求解单个方程？

最佳答案

此答案比较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

如果您比较这些方法，最好听到有关实际问题时机的报告。谢谢。

关于julia - 如何在Julia中快速求解多项式方程式？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/43456595/