Lagrange Interpolation: Difference between revisions
Content added Content deleted
(add RPL + typo in task definition) |
|||
Line 35: | Line 35: | ||
p = Lagrange(1⋅ℓ_0(x) + 4⋅ℓ_1(x) + 1⋅ℓ_2(x) + 5⋅ℓ_3(x)) |
p = Lagrange(1⋅ℓ_0(x) + 4⋅ℓ_1(x) + 1⋅ℓ_2(x) + 5⋅ℓ_3(x)) |
||
Polynomial(p) = Polynomial(-21.0 + 35.83333333333333*x - 16.0*x^2 + 2.1666666666666665*x^3) |
Polynomial(p) = Polynomial(-21.0 + 35.83333333333333*x - 16.0*x^2 + 2.1666666666666665*x^3) |
||
</pre> |
|||
=== Rational base polynomial answer === |
|||
<syntaxhighlight lang="julia">using Polynomials |
|||
function Lagrange(pts::Vector{Vector{Int}}) |
|||
xs = first.(pts) |
|||
cs = last.(pts) |
|||
n = length(xs) |
|||
n == 1 && return Polynomial(cs[1]) |
|||
arr = ones(Int, n) |
|||
for j in 2:n |
|||
for k in 1:j |
|||
arr[k] = (xs[k] - xs[j]) * arr[k] |
|||
end |
|||
arr[j] = prod(xs[j] - xs[i] for i in 1:(j - 1)) |
|||
end |
|||
ws = 1 .// arr |
|||
q = Polynomial(0 // 1) |
|||
x = variable(q) |
|||
for i in eachindex(ws) |
|||
m = prod(x - xs[j] for j in eachindex(xs) if j != i) |
|||
q += m * ws[i] * cs[i] |
|||
end |
|||
return q |
|||
end |
|||
const testpoints = [[1, 1], [2, 4], [3, 1], [4, 5]] |
|||
@show Lagrange(testpoints) |
|||
</syntaxhighlight>{{out}} |
|||
<pre> |
|||
Lagrange(testpoints) = Polynomial(-21//1 + 215//6*x - 16//1*x^2 + 13//6*x^3) |
|||
</pre> |
</pre> |
||
Revision as of 04:50, 8 September 2023
Lagrange Interpolation is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
The task is to implement the Lagrange Interpolation formula and use it to solve the example problem to find a polynomial P of degree<4 satisfying P(1)=1 P(2)=4 P(3)=1 P(4)=5 as described at [1]
- Related task Curve that touches three points
F#
// Lagrange Interpolation. Nigel Galloway: September 5th., 2023
let symbol=MathNet.Symbolics.SymbolicExpression.Variable
let qi=MathNet.Symbolics.SymbolicExpression.FromInt32
let eval (g:MathNet.Symbolics.SymbolicExpression) x=let n=Map["x",MathNet.Symbolics.FloatingPoint.Real x] in MathNet.Symbolics.Evaluate.evaluate n g.Expression
let fN g=let x=symbol "x" in g|>List.fold(fun z c->(x-c)*z)(qi 1)
let fG(n,g)=let n,g=n|>List.map qi,g|>List.map qi in List.map2(fun i g->i,g,n|>List.except [i]) n g
let LIF n=fG n|>List.sumBy(fun(ci,bi,c)->bi*(fN c)/(c|>List.fold(fun z c->(ci-c)*z)(qi 1)))
printfn $"%s{LIF([1;2;3;4],[1;4;1;5]).Expand().ToString()}"
- Output:
-21 + 215/6*x - 16*x^2 + 13/6*x^3
Julia
Note the Polynomials module prints polynomials in order from degree 0 to n rather than n to zero degree.
using Polynomials, SpecialPolynomials
const pts = [[1, 1], [2, 4], [3, 1], [4, 5]]
const xs = first.(pts)
const ys = last.(pts)
const p = Lagrange(xs, ys)
@show p Polynomial(p)
- Output:
p = Lagrange(1⋅ℓ_0(x) + 4⋅ℓ_1(x) + 1⋅ℓ_2(x) + 5⋅ℓ_3(x)) Polynomial(p) = Polynomial(-21.0 + 35.83333333333333*x - 16.0*x^2 + 2.1666666666666665*x^3)
Rational base polynomial answer
using Polynomials
function Lagrange(pts::Vector{Vector{Int}})
xs = first.(pts)
cs = last.(pts)
n = length(xs)
n == 1 && return Polynomial(cs[1])
arr = ones(Int, n)
for j in 2:n
for k in 1:j
arr[k] = (xs[k] - xs[j]) * arr[k]
end
arr[j] = prod(xs[j] - xs[i] for i in 1:(j - 1))
end
ws = 1 .// arr
q = Polynomial(0 // 1)
x = variable(q)
for i in eachindex(ws)
m = prod(x - xs[j] for j in eachindex(xs) if j != i)
q += m * ws[i] * cs[i]
end
return q
end
const testpoints = [[1, 1], [2, 4], [3, 1], [4, 5]]
@show Lagrange(testpoints)
- Output:
Lagrange(testpoints) = Polynomial(-21//1 + 215//6*x - 16//1*x^2 + 13//6*x^3)
RPL
[[1 2 3 4][1 4 1 5]] LAGRANGE
- Output:
1: '(13*X^3-96*X^2+215*X-126)/6'
Wren
A polynomial is represented here by a list of coefficients from the lowest to the highest degree. However, the library methods which deal with polynomials expect the coefficients to be presented from highest to lowest degree so we therefore need to reverse the list before calling these methods.
import "./dynamic" for Tuple
import "./math" for Math
import "./fmt" for Fmt
var Point = Tuple.create("Point", ["x", "y"])
// Add two polynomials.
var add = Fn.new { |p1, p2|
var m = p1.count
var n = p2.count
var sum = List.filled(m.max(n), 0)
for (i in 0...m) sum[i] = p1[i]
for (i in 0...n) sum[i] = sum[i] + p2[i]
return sum
}
// Multiply two polynmials.
var multiply = Fn.new { |p1, p2|
var m = p1.count
var n = p2.count
var prod = List.filled(m + n - 1, 0)
for (i in 0...m) {
for (j in 0...n) prod[i+j] = prod[i+j] + p1[i] * p2[j]
}
return prod
}
// Multiply a polynomial by a scalar.
var scalarMultiply = Fn.new { |poly, x| poly.map { |coef| coef * x }.toList }
// Divide a polynomial by a scalar.
var scalarDivide = Fn.new { |poly, x| scalarMultiply.call(poly, 1/x) }
// Returns the Lagrange interpolating polynomial which passes through a list of points.
var lagrange = Fn.new { |pts|
var c = pts.count
var polys = List.filled(c, null)
for (i in 0...c) {
var poly = [1]
for (j in 0...c) {
if (i == j) continue
poly = multiply.call(poly, [-pts[j].x, 1])
}
var d = Math.evalPoly(poly[-1..0], pts[i].x)
polys[i] = scalarDivide.call(poly, d)
}
var sum = [0]
for (i in 0...c) {
polys[i] = scalarMultiply.call(polys[i], pts[i].y)
sum = add.call(sum, polys[i])
}
return sum
}
var pts = [
Point.new(1, 1),
Point.new(2, 4),
Point.new(3, 1),
Point.new(4, 5)
]
var lip = lagrange.call(pts)
Fmt.pprint("$f", lip[-1..0], "", "x")
- Output:
2.166667x³ - 16.000000x² + 35.833333x - 21.000000