# Interpolation and Differentiation

All algorithms use the baryentric weights as described in the book "Implementing Spectral Methods for PDEs" by David Kopriva.

- Interpolation can be performed either with `interpolate(dest, u, basis)` or via an 
  `interpolation_matrix(dest, basis)`.
- In a nodal basis, the derivative matrix is available as `basis.D`. Alternatively, `derivative_at(x, u, basis)`
  can be used.


```julia
using Revise
using PolynomialBases
using LaTeXStrings, Plots

# define nodal bases
p = 5 # polynomial degree
basis1 = LobattoLegendre(p)
basis2 = GaussLegendre(p)

# the function that will be interpolated
ufunc(x) = sinpi(x); uprim(x) = π*cospi(x)
#ufunc(x) = 1 / (1 + 25x^2); uprim(x) = -ufunc(x)^2*50x

for basis in (basis1, basis2)
    u = ufunc.(basis.nodes)

    xplot = range(-1, stop=1, length=500)
    uplot = interpolate(xplot, u, basis)

    fig1 = plot(xplot, ufunc.(xplot), label="u", xguide=L"x", yguide=L"u")
    plot!(fig1, xplot, uplot, label=L"\mathrm{I}(u)")

    fig2 = plot(xplot, uprim.(xplot), label="u'", xguide=L"x", yguide=L"u'")
    plot!(fig2, xplot, interpolate(xplot, basis.D*u, basis), label=L"\mathrm{I}(u)'")

    display(basis)
    display(plot(fig1, fig2))
end
```

# Integration

The nodes and weights are from [FastGaussQuadrature.jl](https://github.com/ajt60gaibb/FastGaussQuadrature.jl).


```julia
using Revise
using PolynomialBases
using LaTeXStrings, Plots; pyplot()

ufunc(x) = sinpi(x)^6

function compute_error(p, basis_type)
    basis = basis_type(p)
    u = ufunc.(basis.nodes)
    abs(5/8 - integrate(u, basis))
end

ps = 1:23
scatter(ps, compute_error.(ps, LobattoLegendre), label="Lobatto", xguide=L"p", yguide="Error", yaxis=:log10)
scatter!(ps, compute_error.(ps, GaussLegendre), label="Gauss")
```

# Evaluation of Orthogonal Polynomials

## Legendre Polynomials

Legendre poylnomials $P_p$ are evaluated as `legendre(x, p)` using the three term recursion formula.


```julia
using Revise
using PolynomialBases
using LaTeXStrings, Plots; pyplot()

x = range(-1, stop=1, length=10^3)
fig = plot(xguide=L"x")
for p in 0:5
    plot!(fig, x, legendre.(x, p), label="\$ P_$p \$")
end
fig
```

## Gegenbauer Polynomials

Gegenbauer poylnomials $C_p^{(\alpha)}$ are evaluated as `gegenbauer(x, p, α)` using the three term recursion formula.


```julia
using Revise
using PolynomialBases
using LaTeXStrings, Plots; pyplot()

α = 0.5
x = range(-1, stop=1, length=10^3)
fig = plot(xguide=L"x")
for p in 0:5
    plot!(fig, x, gegenbauer.(x, p, α), label="\$ C_$p^{($α)} \$")
end
fig
```

## Jacobi Polynomials

Jacobi poylnomials $P_p^{\alpha,\beta}$ are evaluated as `jacobi(x, p, α, β)` using the three term recursion formula.


```julia
using Revise
using PolynomialBases
using LaTeXStrings, Plots; pyplot()

α, β = -0.5, -0.5
#α, β = 0, 0
x = range(-1, stop=1, length=10^3)
fig = plot(xguide=L"x")
for p in 0:5
    plot!(fig, x, jacobi.(x, p, α, β), label="\$ P_$p^{$α, $β} \$")
end
fig
```

## Hermite Polynomials

Hermite poylnomials $H_p$ are evaluated as `hermite(x, p)` using the three term recursion formula.


```julia
using Revise
using PolynomialBases
using LaTeXStrings, Plots; pyplot()

x = range(-2, stop=2, length=10^3)
fig = plot(xguide=L"x")
for p in 0:4
    plot!(fig, x, hermite.(x, p), label="\$ H_$p \$")
end
fig
```

## Hahn Polynomials


```julia
using Revise
using PolynomialBases
using LaTeXStrings, Plots; pyplot()

α, β, N = -0.5, -0.5, 20
x = range(0, stop=N, length=10^3)
fig = plot(xguide=L"x")
for p in 0:5
    plot!(fig, x, hahn.(x, p, α, β, N), label="\$ Q_$p(x; $α, $β, $N) \$")
end
fig
```

# Symbolic Computations

Symbolic computations using [SymPy.jl](https://github.com/JuliaPy/SymPy.jl) and [SymEngine.jl](https://github.com/symengine/symengine) are supported.


```julia
using Revise
using PolynomialBases
import SymPy
import SymEngine

x_sympy = SymPy.symbols("x")
x_symengine = SymEngine.symbols("x")

legendre(x_sympy, 6) |> SymPy.expand |> display
legendre(x_symengine, 6) |>  SymEngine.expand |> display

GaussLegendre(2, SymPy.Sym).D |> display
GaussLegendre(2, SymEngine.Basic).D |> display
```

# Modal Matrices

The Vandermonde matrix $V$ can be computed as `legendre_vandermonde(basis)` and used to transform coeficients in a modal basis of Legendre polynomials to a nodal basis.


```julia
using Revise
using PolynomialBases
using LinearAlgebra

p = 7
basis = GaussLegendre(p)
V = legendre_vandermonde(basis)

# the modal derivative matrix
Dhat = legendre_D(p)
# they should be equal
norm( basis.D - V * Dhat / V)
```

Similarly, the Vandermonde matrix with respect to the Jacobi polynomials is given as `jacobi_vandermonde(basis, α, β)`.


```julia
using Revise
using PolynomialBases

p = 3
α, β = -0.5, -0.5
basis = GaussJacobi(p, α, β)
V = jacobi_vandermonde(basis, α, β)
```

a64a97b49a9b9005f76f93c9bda60dc86e9e4bec
ipynb
Jupyter Notebook
notebooks/Tutorial.ipynb
ranocha/PolynomialBases.jl
845531c8eee63f2c1608badd980872796fabbfc9
MIT
2019-09-11T18:12:07.000Z
2022-03-13T21:08:56.000Z
notebooks/Tutorial.ipynb
ranocha/PolynomialBases.jl
845531c8eee63f2c1608badd980872796fabbfc9
MIT
2018-02-15T06:32:55.000Z
2022-01-11T15:59:07.000Z
notebooks/Tutorial.ipynb
ranocha/PolynomialBases.jl
845531c8eee63f2c1608badd980872796fabbfc9
MIT
2018-02-26T18:34:02.000Z
2020-02-08T11:01:42.000Z
Qwen/Qwen-72B
1. YES
2. YES