How to fit a polynomial with a pre-determined degree to a set of points?

General Usage
, , , ,
3

The thing is that ApproxFun needs to know the values of f at specific values of x. For a Chebyshev approximation, if you happen to know the function values at the Chebyshev nodes there’s no problem:

using CairoMakie
using ApproxFun

f(x) = sin(2*pi*x)

# Get 30 Chebyshev nodes in the interval [0, 5]
# (these nodes are not evenly spaced)
S = Chebyshev(0..5)
x = points(S, 30)

# Suppose we know the value of f at these points
y = f.(x)

# We can calculate the corresponding Chebyshev coefficients
c = transform(S, y)
# And build the approximation of f
fa = Fun(S, c)

fig = plot(x, y, figure=(; resolution=(500, 200)))
plot!(domain(S), fa)
fig

image

More probably, you have data points on an even grid or a random grid. You can still find the coefficients of the Chebyshev polynomials that go through these points using linear regression as shown in the FAQ:

# Points on a regular grid
x = range(0, 5, length=30)
y = f.(x)

# Function space for approximations on the given interval
S = Chebyshev(x[begin]..x[end])

# Build design matrix for linear regression by evaluating each
# polynomial at the given x values
X = hcat((Fun(S, [zeros(k-1); 1]).(x) for k in 1:length(x))...)

# Use linear regression to find the Chebyshev coefficients
c = X\y
fa = Fun(S, c)

fig = plot(x, y, figure=(; resolution=(500, 200)))
plot!(domain(S), fa)
fig

image

Though according to the FAQ for numberical stability it’s better to use an overdetermined system (i.e. replace length(x) by a smaller number in the code above).

Another option proposed here is to make a first approximation with the AAA algorithm and use the result (a rational function) to calculate values at the Chebyshev nodes:

using BaryRational

# Using >30 points to avoid Froissart doublets (poles near zeros)
x = range(0, 5, length=50)
y = f.(x)

# Rational approximation
f_aaa = aaa(x, y)

# Approximations of f at the Chebyshev nodes (we choose n=100)
x_cheb = points(S, 100)
y_cheb = f_aaa.(x_cheb)

# Chebyshev approximation of f
c = transform(S, y_cheb)
fa = Fun(S, c)

fig = plot(x, y, figure=(; resolution=(500, 200)))
plot!(0..5, fa)
fig

image

3 Likes