In [1]:
using Plots, LinearAlgebra

Vandermonde Matrix¶

The Vandermonde matrix naturally arises in approximating a function using a polynomial:

$$ f(x) = a_0 + a_1 x + a_2 x^2 + \ldots + a_{N-1} x^{N-1}. $$

By sampling the functions at $N$ distinct points, the coefficients $a_n$'s satisfy

$$ \begin{aligned} V \begin{bmatrix} a_0 \\ a_1 \\ \ldots \\ a_{N-1} \end{bmatrix} = \begin{bmatrix} f(x_0) \\ f(x_1) \\ \ldots \\ f(x_{N-1}) \end{bmatrix} \end{aligned}, $$ where

$$ \begin{aligned} V = \begin{bmatrix} 1 & x_0 & x_0^2 & \ldots & x_0^{N-1} \\ 1 & x_1 & x_1^2 & \ldots & x_1^{N-1} \\ 1 & x_2 & x_2^2 & \ldots & x_2^{N-1} \\ & & \ldots & & \\ 1 & x_{N-1} & x_{N-1}^2 & \ldots & x_{N-1}^{N-1} \end{bmatrix} \end{aligned} $$

In [2]:
function Vandermonde(xv)
    
    ndim = length(xv)
    vmat = zeros(ndim, ndim)
    
    for i1 in 1:ndim, i2 in 1:ndim
			vmat[i1,i2] = xv[i1]^(i2-1)
	end
    return vmat
    
end

function Vinterp(f, xval; a=-1, b=1, ndim=25)

    xv = range(a, b, length=ndim)
    coeffs = Vandermonde(xv) \ f.(xv)
    ret = sum( [coeffs[j] * xval^(j-1) for j in 1:ndim] )
    
    return ret
    
end
Out[2]:
Vinterp (generic function with 1 method)
In [3]:
f = x -> 1/(1+x^2)

_f1 = x -> Vinterp(f, x, ndim=5)
xv = range(-1.2, 1.2, length=20)

p1 = plot(xv, f, label="exact")
scatter!(xv, _f1, label="interpolation (5th)")
Out[3]:

Quadrature Formula¶

$$ \int_{-1}^{1} dx f(x) \approx \sum_{n={\rm even}} \frac{2}{n+1} a_n = \sum_{m} y_m \sum_{n={\rm even}} \frac{2}{n+1} V^{-1}_{nm} $$

In [4]:
function quad_coeff(i1; norder=3)

    vv = Vandermonde(range(-1, 1, length=norder))
    
    ret = 0.0

    for j1 in 1:norder
        if j1%2 == 0
            continue
        else
            ret += inv(vv)[j1, i1] * (2//(j1))
        end
    end

    return ret
end
Out[4]:
quad_coeff (generic function with 1 method)
In [5]:
for nn in 2:5
    println(quad_coeff.(1:nn, norder=nn))
end
[1.0, 1.0]
[0.3333333333333333, 1.3333333333333335, 0.3333333333333333]
[0.25, 0.7499999999999999, 0.7500000000000002, 0.24999999999999992]
[0.15555555555555547, 0.711111111111111, 0.26666666666666705, 0.711111111111111, 0.1555555555555556]