Общая функция для решения полиномиальных корней n-го порядка в Юлии

Все,

Я только начал играть с языком Джулия и мне это очень нравится. В конце 3-го урока есть интересная проблема: обобщить квадратную формулу так, чтобы она решала корни любого полиномиального уравнения n-го порядка.

Это показалось мне (а) интересной проблемой программирования и (б) интересной проблемой Юлии. Кто-нибудь там решил это? Для справки вот код Джулии с парой игрушечных примеров. Опять же, идея состоит в том, чтобы сделать это универсальным для любого полинома 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)
Источник

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 также имеет поддержку для вычисления произвольной точности. Это полезно для решения уравнения, которое не может быть решено с двойной точностью. Например

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)
Источник

Алгебраических формул для общих многочленов пятой степени и выше не существует ( см. Здесь). Теоретически, вы могли бы продолжить использовать ту же методологию для решений кубических и квартичных, но даже это было бы большой тяжелой работой, учитывая очень громоздкие формулы для корней квартичных. Вы также можете использовать CAS, такой как SymPy, чтобы узнать эти формулы.

Источник
Другие вопросы по тегам