copy over from OrthogonalPolynomialsQuasi.jl

Author: dlfivefifty

.github/workflows/CompatHelper.yml

@@ -0,0 +1,25 @@
+name: CompatHelper
+
+on:
+  schedule:
+    - cron: '00 * * * *'
+
+jobs:
+  build:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        julia-version: [1.2.0]
+        julia-arch: [x86]
+        os: [ubuntu-latest]
+    steps:
+      - uses: julia-actions/setup-julia@latest
+        with:
+          version: ${{ matrix.julia-version }}
+      - name: Install dependencies
+        run: julia -e 'using Pkg; Pkg.add(Pkg.PackageSpec(name = "CompatHelper", url = "https://github.com/bcbi/CompatHelper.jl.git"))'
+      - name: CompatHelper.main
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          JULIA_DEBUG: CompatHelper
+        run: julia -e 'using CompatHelper; CompatHelper.main()'

.github/workflows/TagBot.yml

@@ -0,0 +1,17 @@
+name: TagBot
+on:
+  issue_comment:  # THIS BIT IS NEW
+    types:
+      - created
+  workflow_dispatch:
+jobs:
+  TagBot:
+    # THIS 'if' LINE IS NEW
+    if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot'
+    # NOTHING BELOW HAS CHANGED
+    runs-on: ubuntu-latest
+    steps:
+      - uses: JuliaRegistries/TagBot@v1
+        with:
+          token: ${{ secrets.GITHUB_TOKEN }}
+          ssh: ${{ secrets.DOCUMENTER_KEY }}

.github/workflows/ci.yml

@@ -0,0 +1,42 @@
+name: CI
+on:
+  - push
+  - pull_request
+jobs:
+  test:
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version:
+          - '1.5'
+          - 'nightly'
+        os:
+          - ubuntu-latest
+          - macOS-latest
+          - windows-latest
+        arch:
+          - x64
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: actions/cache@v1
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v1
+        with:
+          file: lcov.info

.gitignore

@@ -22,3 +22,4 @@ docs/site/
 # committed for packages, but should be committed for applications that require a static
 # environment.
 Manifest.toml
+.DS_Store

Project.toml

@@ -0,0 +1,52 @@
+name = "ClassicalOrthogonalPolynomials"
+uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
+authors = ["Sheehan Olver <[email protected]>"]
+version = "0.1.0"
+
+[deps]
+ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
+BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
+BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
+ContinuumArrays = "7ae1f121-cc2c-504b-ac30-9b923412ae5c"
+DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
+FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c"
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
+InfiniteLinearAlgebra = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
+IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
+LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+
+[compat]
+ArrayLayouts = "0.5.1"
+BandedMatrices = "0.16"
+BlockArrays = "0.14"
+BlockBandedMatrices = "0.10"
+ContinuumArrays = "0.5"
+DomainSets = "0.4, 0.5"
+FFTW = "1.1"
+FastGaussQuadrature = "0.4.3"
+FastTransforms = "0.11"
+FillArrays = "0.11"
+InfiniteArrays = "0.9"
+InfiniteLinearAlgebra = "0.4.6"
+IntervalSets = "0.3.1, 0.4, 0.5"
+LazyArrays = "0.20"
+QuasiArrays = "0.4"
+SpecialFunctions = "0.10, 1"
+julia = "1.5"
+
+[extras]
+Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
+SemiseparableMatrices = "f8ebbe35-cbfb-4060-bf7f-b10e4670cf57"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Base64", "Test", "ForwardDiff", "SemiseparableMatrices", "LazyBandedMatrices"]

README.md

@@ -1,2 +1,89 @@
 # ClassicalOrthogonalPolynomials.jl
 A Julia package for classical orthogonal polynomials and expansions
+
+[![Build Status](https://github.com/JuliaApproximation/ClassicalOrthogonalPolynomials.jl/workflows/CI/badge.svg)](https://github.com/JuliaApproximation/FastGaussQuadrature.jl/actions)
+[![codecov](https://codecov.io/gh/JuliaApproximation/ClassicalOrthogonalPolynomials.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/JuliaApproximation/ClassicalOrthogonalPolynomials.jl)
+
+
+This package implements classical orthogonal polynomials as quasi-arrays where one one axes is continuous and the other axis is discrete (countably infinite), as implemented in [QuasiArrays.jl](https://github.com/JuliaApproximation/QuasiArrays.jl) and  [ContinuumArrays.jl](https://github.com/JuliaApproximation/ContinuumArrays.jl).  
+```julia
+julia> using ClassicalOrthogonalPolynomials, ContinuumArrays
+
+julia> P = Legendre(); # Legendre polynomials
+
+julia> size(P) # uncountable ∞ x countable ∞
+(ℵ₁, ∞)
+
+julia> axes(P) # essentially (-1..1, 1:∞), Inclusion plays the same role as Slice
+(Inclusion(-1.0..1.0 (Chebyshev)), OneToInf())
+
+julia> P[0.1,1:10] # [P_0(0.1), …, P_9(0.1)]
+10-element Array{Float64,1}:
+  1.0                
+  0.1                
+ -0.485              
+ -0.14750000000000002
+  0.3379375          
+  0.17882875         
+ -0.2488293125       
+ -0.19949294375000004
+  0.180320721484375  
+  0.21138764183593753
+
+julia> @time P[range(-1,1; length=10_000), 1:10_000]; # construct 10_000^2 Vandermonde matrix
+  1.624796 seconds (10.02 k allocations: 1.491 GiB, 6.81% gc time)
+```
+This also works for associated Legendre polynomials as weighted Ultraspherical polynomials:
+```julia
+julia> associatedlegendre(m) = ((-1)^m*prod(1:2:(2m-1)))*(UltrasphericalWeight((m+1)/2).*Ultraspherical(m+1/2))
+associatedlegendre (generic function with 1 method)
+
+julia> associatedlegendre(2)[0.1,1:10]
+10-element Array{Float64,1}:
+   2.9699999999999998
+   1.4849999999999999
+  -6.9052500000000006
+  -5.041575          
+  10.697754375       
+  10.8479361375      
+ -13.334647528125    
+ -18.735466024687497 
+  13.885467170308594 
+  28.220563705988674 
+```
+
+## p-Finite Element Method
+
+The language of quasi-arrays gives a natural framework for constructing p-finite element methods. The convention
+is that adjoint-products are understood as inner products over the axes with uniform weight. Thus to solve Poisson's equation
+using its weak formulation with Dirichlet conditions we can expand in a weighted Jacobi basis:
+```julia
+julia> P¹¹ = Jacobi(1.0,1.0); # Quasi-matrix of Jacobi polynomials
+
+julia> w = JacobiWeight(1.0,1.0); # quasi-vector correspoinding to (1-x^2)
+
+julia> w[0.1] ≈ (1-0.1^2)
+true
+
+julia> S = w .* P¹¹; # Quasi-matrix of weighted Jacobi polynomials
+
+julia> D = Derivative(axes(S,1)); # quasi-matrix corresponding to derivative
+
+julia> Δ = (D*S)'*(D*S) # weak laplacian corresponding to inner products of weighted Jacobi polynomials
+∞×∞ LazyArrays.ApplyArray{Float64,2,typeof(*),Tuple{Adjoint{Int64,BandedMatrices.BandedMatrix{Int64,Adjoint{Int64,InfiniteArrays.InfStepRange{Int64,Int64}},InfiniteArrays.OneToInf{Int64}}},LazyArrays.BroadcastArray{Float64,2,typeof(*),Tuple{LazyArrays.BroadcastArray{Float64,1,typeof(/),Tuple{Int64,InfiniteArrays.InfStepRange{Int64,Int64}}},BandedMatrices.BandedMatrix{Int64,Adjoint{Int64,InfiniteArrays.InfStepRange{Int64,Int64}},InfiniteArrays.OneToInf{Int64}}}}}} with indices OneToInf()×OneToInf():
+ 2.66667   ⋅     ⋅        ⋅        ⋅        ⋅        ⋅        ⋅      …  
+  ⋅       6.4    ⋅        ⋅        ⋅        ⋅        ⋅        ⋅         
+  ⋅        ⋅   10.2857    ⋅        ⋅        ⋅        ⋅        ⋅         
+  ⋅        ⋅     ⋅      14.2222    ⋅        ⋅        ⋅        ⋅         
+  ⋅        ⋅     ⋅        ⋅      18.1818    ⋅        ⋅        ⋅         
+  ⋅        ⋅     ⋅        ⋅        ⋅      22.1538    ⋅        ⋅      …  
+  ⋅        ⋅     ⋅        ⋅        ⋅        ⋅      26.1333    ⋅         
+  ⋅        ⋅     ⋅        ⋅        ⋅        ⋅        ⋅      30.1176     
+  ⋅        ⋅     ⋅        ⋅        ⋅        ⋅        ⋅        ⋅         
+  ⋅        ⋅     ⋅        ⋅        ⋅        ⋅        ⋅        ⋅         
+  ⋅        ⋅     ⋅        ⋅        ⋅        ⋅        ⋅        ⋅      …  
+  ⋅        ⋅     ⋅        ⋅        ⋅        ⋅        ⋅        ⋅         
+ ⋮                                         ⋮                         ⋱  
+```
+
+

examples/beam.jl

@@ -0,0 +1,66 @@
+using ClassicalOrthogonalPolynomials, ContinuumArrays, DifferentialEquations, Plots
+
+###
+# Heat
+###
+
+S = JacobiWeight(1.0,1.0) .* Jacobi(1.0,1.0)
+D = Derivative(axes(S,1))
+Δ = -((D*S)'*(D*S))
+M = S'S
+
+function evolution(u, (M,A), t) 
+    M\(A*u)
+end
+n = 50
+u0 = [[1,2,3]; zeros(n-3)]
+prob = ODEProblem(evolution,u0,(0.0,3.0),(cholesky(Symmetric(M[1:n,1:n])),Δ[1:n,1:n]))
+@time u = solve(prob,TRBDF2()); u(1.0)
+
+###
+# Heat natural
+###
+L = LinearSpline(range(-1,1;length=2))
+P = apply(hcat,L,S)
+ApplyArray(hcat,Legendre() \ (D*L), Legendre() \ (D*S))
+
+(D*L)'*(D*L)
+
+(D*S)'*(D*S)
+
+((D*S)'*(D*L)).args[3]
+(D*L)'*(D*S)
+H = (D*L).args[1]
+P = Legendre()
+(P\H)' * P'P * (P \ (D*S))
+
+
+
+###
+# Beam
+###
+
+# [u, u_t]_t = [0 I; M] * [u, u_t]
+
+S = JacobiWeight(2.0,2.0) .* Jacobi(2.0,2.0)
+D = Derivative(axes(S,1))
+Δ² = (D^2*S)'*(D^2*S)
+M = S'S
+function wave(uv, (M,A), t) 
+    n = size(M,1)
+    u,v = uv[1:n],uv[n+1:end]
+    [v; M\(A*u)]
+end
+n = 30
+u0 = [1; 2; 3; zeros(n-3)]
+prob = ODEProblem(wave,[u0; zeros(n)],(0.0,3.0),(cholesky(Symmetric(M[1:n,1:n])),-(Δ²)[1:n,1:n]))
+@time u = solve(prob,TRBDF2()); u(1.0)
+
+g = range(-1,1;length=200)
+V = S[g,1:n]
+@gif for t in 0.0:0.01:3.0
+    plot(g, V*u(t)[1:n]; ylims=(-5,5))
+end
+
+
+prob = ODEProblem(wave,[u0; zeros(n)],(0.0,3.0),(qr(M),-(Δ²)))

examples/collocation.jl

@@ -0,0 +1,21 @@
+using ContinuumArrays, QuasiArrays, FillArrays, InfiniteArrays
+import QuasiArrays: Inclusion, QuasiDiagonal
+
+using Plots; pyplot();
+
+S = LinearSpline(range(0,1; length=10))
+xx = range(0,1; length=1000)
+
+S = Jacobi(true,true)
+
+
+P = Legendre()
+X = QuasiDiagonal(Inclusion(-1..1))
+
+@test X[-1:0.1:1,-1:0.1:1] == Diagonal(-1:0.1:1)
+
+axes(X)
+J = pinv(P)*X*P
+
+J - I
+Vcat(Hcat(1, Zeros(1,∞)), J)

examples/golubwelsch.jl

@@ -0,0 +1,5 @@
+using ClassicalOrthogonalPolynomials, FastGaussQuadrature
+
+P = Legendre()
+x = axes(P,1)
+x .* P

examples/ultrasphericalspectralmethod.jl

@@ -0,0 +1,14 @@
+using ContinuumArrays, FillArrays, InfiniteArrays, Plots
+
+T = Chebyshev()
+C = Ultraspherical(2)
+D = Derivative(axes(T,1))
+A = (C\(D^2*T))+100(C\T)
+
+n = 100
+c = [T[[-1,1],1:n]; A[1:n-2,1:n]] \ [1;2;zeros(n-2)]
+u = T*Vcat(c,Zeros(∞))
+
+xx = range(-1,1;length=1000)
+plot(xx,u[xx])
+