全部,
我刚刚开始玩转Julia语言,并且非常喜欢它。在第3篇教程的最后,有一个有趣的问题:对二次公式进行通用化,使其可求解任何n-order polynomial equation的根。
这让我感到震惊(a)一个有趣的编程问题,(b)一个有趣的Julia问题。有没有人解决了这个问题?作为引用,这是带有两个玩具示例的Julia代码。再次,该想法是使该泛型适用于任何n阶多项式。
干杯,
亚伦
function derivative(f)
return function(x)
# pick a small value for h
h = x == 0 ? sqrt(eps(Float64)) : sqrt(eps(Float64)) * x
# floating point arithmetic gymnastics
xph = x + h
dx = xph - x
# evaluate f at x + h
f1 = f(xph)
# evaluate f at x
f0 = f(x)
# divide the difference by h
return (f1 - f0) / dx
end
end
function quadratic(f)
f1 = derivative(f)
c = f(0.0)
b = f1(0.0)
a = f(1.0) - b - c
return (-b + sqrt(b^2 - 4a*c + 0im))/2a, (-b - sqrt(b^2 - 4a*c + 0im))/2a
end
quadratic((x) -> x^2 - x - 2)
quadratic((x) -> x^2 + 2)
最佳答案
软件包 PolynomialRoots.jl
提供函数roots()
来查找任何阶数的多项式的所有(实数和复数)根。唯一的强制性参数是多项式系数按升序排列的数组。
例如,为了找到根源
6x^5 + 5x^4 + 3x^2 + 2x + 1
加载包(
using PolynomialRoots
)后,您可以使用julia> roots([1, 2, 3, 4, 5, 6])
5-element Array{Complex{Float64},1}:
0.294195-0.668367im
-0.670332+2.77556e-17im
0.294195+0.668367im
-0.375695-0.570175im
-0.375695+0.570175im
该软件包是本文描述的寻根算法的Julia实现:http://arxiv.org/abs/1203.1034
PolynomialRoots.jl
还支持任意精度计算。这对于求解无法以 double 求解的方程式很有用。例如julia> r = roots([94906268.375, -189812534, 94906265.625]);
julia> (r[1], r[2])
(1.0000000144879793 - 0.0im,1.0000000144879788 + 0.0im)
给出错误的多项式结果,而是以任意精度传递输入数组会强制提供正确答案的任意精度计算(请参阅https://en.wikipedia.org/wiki/Loss_of_significance):
julia> r = roots([BigFloat(94906268.375), BigFloat(-189812534), BigFloat(94906265.625)]);
julia> (Float64(r[1]), Float64(r[2]))
(1.0000000289759583,1.0)