CODEWITHSUNDEEP
julianumerical-methods
jakekap7
jakekap7

Reputation: 35

Algorithm for Horner's Method with base points in Julia

I'm attempting to write an algorithm in Julia for Horner's method with base points (b1, b2, etc.) to evaluate the polynomial: y = p(x) = c1 + c2(x-b1) + c3(x-b1)(x-b2) + c4(x-b1)(x-b2)(x-b3), where the coefficients are represented by "c" and the base points are represented by "b". The polynomial that I need to evaluate is: p(x) = 10 - (21/2)(x+2) + (45/6)(x+2)(x) - (7/3)(x+2)(x)(x-1).

Here is my code so far:

function hornerpolybase(x,c,b)
    
    # assign y to the first coefficient
    # make y the same type as input x
    y = c[1] * one(x) 
    
    n = length(c)
    
    for i = 2:n
        y = y + c[i] * (x - b[i-1])   
    end
    
    return y
end

My issue is that after the second iteration of the for loop, y = c1 + c2(x-b1) + c3(x-b2), whereas the equation should be y = c1 + c2(x-b1) + c3 * (x-b1)(x-b2). The problem repeats in the third iteration. I'm thinking that possibly an "if, then" conditional statement within the for loop should solve this issue, but I'm completely stuck. I would appreciate any help.

Upvotes: 0

Views: 265

Answers (2)

user1196549
user1196549

Reputation:

The standard Horner's scheme goes

C0 + C1 X + C2 X² + C3 X³ = C0 + X (C1 + C2 X + C3 X²)
                          = C0 + X (C1 + X (C2 + C3 X))

So what you compute is, in this order

C2 + X C3
C1 + X (C2 + C3 X)
C0 + X (C1 + X (C2 + C3 X))

Using an intermediate variable,

Y = C3
Y = C2 + X Y
Y = C1 + X Y
Y = C0 + X Y

In your case, there are also offsets:

Y = C3
Y = C2 + (X - X2) Y
Y = C1 + (X - X1) Y
Y = C0 + (X - X0) Y

Upvotes: 0

August
August

Reputation: 12558

  • The code isn't fully type stable, despite the way y is initialized. E.g. try @code_warntype hornerpolybase(1, [1,1], [1.0]). You could fix it by including the type of b:

    # assign y to the first coefficient
T = promote_type(typeof(x), eltype(c), eltype(b))
y::T = c[1] # implicit conversion to T

    (evalpoly in Julia actually has the same problem...)

  • You're accessing c in reverse (if c[1] is supposed to correspond with c_1). You need to start with the last coefficient in Horner's method.

  • In each iteration it's y that should be multiplied, not the coefficient:

    y::T = c[end]

n = length(c)
for i = n-1:-1:1
    y = c[i] + y * (x - b[i])   
end

Upvotes: 3

Related Questions