From 5b24742f716e831b0d289eecbdf19421ab3ee9d1 Mon Sep 17 00:00:00 2001 From: Shuhei Ohno Date: Fri, 12 Jul 2024 03:55:43 +0900 Subject: [PATCH] v0.0.1 --- .github/dependabot.yml | 7 + .github/workflows/CI.yml | 71 +++++++ .github/workflows/CompatHelper.yml | 16 ++ .github/workflows/TagBot.yml | 31 +++ .gitignore | 27 +-- Project.toml | 16 ++ README.md | 68 ++++++- dev/deploy.jl | 6 + dev/docs.jl | 3 + dev/pkg.jl | 29 +++ dev/revice.jl | 10 + dev/test.jl | 9 + docs/Project.toml | 3 + docs/make.jl | 24 +++ docs/src/API.md | 12 ++ docs/src/index.md | 72 +++++++ src/FiniteDifferenceMatrices.jl | 302 +++++++++++++++++++++++++++++ test/runtests.jl | 27 +++ 18 files changed, 708 insertions(+), 25 deletions(-) create mode 100644 .github/dependabot.yml create mode 100644 .github/workflows/CI.yml create mode 100644 .github/workflows/CompatHelper.yml create mode 100644 .github/workflows/TagBot.yml create mode 100644 Project.toml create mode 100644 dev/deploy.jl create mode 100644 dev/docs.jl create mode 100644 dev/pkg.jl create mode 100644 dev/revice.jl create mode 100644 dev/test.jl create mode 100644 docs/Project.toml create mode 100644 docs/make.jl create mode 100644 docs/src/API.md create mode 100644 docs/src/index.md create mode 100644 src/FiniteDifferenceMatrices.jl create mode 100644 test/runtests.jl diff --git a/.github/dependabot.yml b/.github/dependabot.yml new file mode 100644 index 0000000..700707c --- /dev/null +++ b/.github/dependabot.yml @@ -0,0 +1,7 @@ +# https://docs.github.com/github/administering-a-repository/configuration-options-for-dependency-updates +version: 2 +updates: + - package-ecosystem: "github-actions" + directory: "/" # Location of package manifests + schedule: + interval: "weekly" diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml new file mode 100644 index 0000000..9324055 --- /dev/null +++ b/.github/workflows/CI.yml @@ -0,0 +1,71 @@ +name: CI +on: + push: + branches: + - main + tags: ['*'] + pull_request: + workflow_dispatch: +concurrency: + # Skip intermediate builds: always. + # Cancel intermediate builds: only if it is a pull request build. + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} +jobs: + test: + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + runs-on: ${{ matrix.os }} + timeout-minutes: 60 + permissions: # needed to allow julia-actions/cache to proactively delete old caches that it has created + actions: write + contents: read + strategy: + fail-fast: false + matrix: + version: + - '1.0' + - '1.7' + os: + - ubuntu-latest + arch: + - x64 + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v2 + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-runtest@v1 + docs: + name: Documentation + runs-on: ubuntu-latest + permissions: + actions: write # needed to allow julia-actions/cache to proactively delete old caches that it has created + contents: write + statuses: write + steps: + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 + with: + version: '1' + - uses: julia-actions/cache@v2 + - name: Configure doc environment + shell: julia --project=docs --color=yes {0} + run: | + using Pkg + Pkg.develop(PackageSpec(path=pwd())) + Pkg.instantiate() + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-docdeploy@v1 + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} + - name: Run doctests + shell: julia --project=docs --color=yes {0} + run: | + using Documenter: DocMeta, doctest + using FiniteDifferenceMatrices + DocMeta.setdocmeta!(FiniteDifferenceMatrices, :DocTestSetup, :(using FiniteDifferenceMatrices); recursive=true) + doctest(FiniteDifferenceMatrices) diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml new file mode 100644 index 0000000..cba9134 --- /dev/null +++ b/.github/workflows/CompatHelper.yml @@ -0,0 +1,16 @@ +name: CompatHelper +on: + schedule: + - cron: 0 0 * * * + workflow_dispatch: +jobs: + CompatHelper: + runs-on: ubuntu-latest + steps: + - name: Pkg.add("CompatHelper") + run: julia -e 'using Pkg; Pkg.add("CompatHelper")' + - name: CompatHelper.main() + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }} + run: julia -e 'using CompatHelper; CompatHelper.main()' diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml new file mode 100644 index 0000000..0cd3114 --- /dev/null +++ b/.github/workflows/TagBot.yml @@ -0,0 +1,31 @@ +name: TagBot +on: + issue_comment: + types: + - created + workflow_dispatch: + inputs: + lookback: + default: "3" +permissions: + actions: read + checks: read + contents: write + deployments: read + issues: read + discussions: read + packages: read + pages: read + pull-requests: read + repository-projects: read + security-events: read + statuses: read +jobs: + TagBot: + if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot' + runs-on: ubuntu-latest + steps: + - uses: JuliaRegistries/TagBot@v1 + with: + token: ${{ secrets.GITHUB_TOKEN }} + ssh: ${{ secrets.DOCUMENTER_KEY }} diff --git a/.gitignore b/.gitignore index 29126e4..5a16984 100644 --- a/.gitignore +++ b/.gitignore @@ -1,24 +1,3 @@ -# Files generated by invoking Julia with --code-coverage -*.jl.cov -*.jl.*.cov - -# Files generated by invoking Julia with --track-allocation -*.jl.mem - -# System-specific files and directories generated by the BinaryProvider and BinDeps packages -# They contain absolute paths specific to the host computer, and so should not be committed -deps/deps.jl -deps/build.log -deps/downloads/ -deps/usr/ -deps/src/ - -# Build artifacts for creating documentation generated by the Documenter package -docs/build/ -docs/site/ - -# File generated by Pkg, the package manager, based on a corresponding Project.toml -# It records a fixed state of all packages used by the project. As such, it should not be -# committed for packages, but should be committed for applications that require a static -# environment. -Manifest.toml +/Manifest.toml +/docs/Manifest.toml +/docs/build/ diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..8ac7183 --- /dev/null +++ b/Project.toml @@ -0,0 +1,16 @@ +name = "FiniteDifferenceMatrices" +uuid = "a7a66f33-e7b8-47af-b618-f9b5bea05f3d" +authors = ["Shuhei Ohno"] +version = "0.0.1" + +[deps] +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" + +[compat] +julia = "1" + +[extras] +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test"] diff --git a/README.md b/README.md index bf71c04..fee321d 100644 --- a/README.md +++ b/README.md @@ -1 +1,67 @@ -# FiniteDifferenceMatrices.jl \ No newline at end of file +```@meta +CurrentModule = FiniteDifferenceMatrices +``` + +# FiniteDifferenceMatrices.jl + +[![Build Status](https://github.com/ohno/FiniteDifferenceMatrices.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/ohno/FiniteDifferenceMatrices.jl/actions/workflows/CI.yml?query=branch%3Amain) +[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://ohno.github.io/FiniteDifferenceMatrices.jl/stable/) +[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://ohno.github.io/FiniteDifferenceMatrices.jl/dev/) + +A Julia package for discrete approximations of differential operators + +## Install + +Run the following code on the to install this package. +```julia +import Pkg; Pkg.add(url="https://github.com/ohno/FiniteDifferenceMatrices.jl.git") +``` + +## Usage + +Run the following code before each use. +```@example index +using FiniteDifferenceMatrices +``` + +A central finite difference of the second-order derivative is +```math +\frac{\mathrm{d}^{2}f(x)}{\mathrm{d} x^{2}} = \frac{f(x+\Delta x) - 2f(x) + f(x-\Delta x)}{\Delta x^{2}} + O(\Delta x^{2}). +``` +```@repl index +fdcoefficient(n=2, m=2, d=:c) +``` + +A discrete approximation of the second-order differential operator is +```math +\frac{\mathrm{d}^2}{\mathrm{d} x^2} +\simeq +\frac{1}{\Delta x^2} +\left(\begin{array}{ccccccc} + -2 & 1 & 0 & \ldots & 0 & 0 & 0 \\ + 1 & -2 & 1 & \ldots & 0 & 0 & 0 \\ + 0 & 1 & -2 & \ldots & 0 & 0 & 0 \\ + \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ + 0 & 0 & 0 & \ldots & -2 & 1 & 0 \\ + 0 & 0 & 0 & \ldots & 1 & -2 & 1 \\ + 0 & 0 & 0 & \ldots & 0 & 1 & -2 +\end{array}\right). +``` +```@repl index +fdmatrix(5, n=2, m=2, d=:c, h=1//1) +``` + +Please see the [API reference](https://ohno.github.io/FiniteDifferenceMatrices.jl/dev/API/) for details and more examples. + +## Developer's Guide + +There are several tools for developers. + +```sh +git clone https://github.com/ohno/FiniteDifferenceMatrices.jl.git +cd FiniteDifferenceMatrices.jl +julia +julia> include("dev/revice.jl") +julia> include("dev/test.jl") +julia> include("dev/docs.jl") +``` diff --git a/dev/deploy.jl b/dev/deploy.jl new file mode 100644 index 0000000..9f791f0 --- /dev/null +++ b/dev/deploy.jl @@ -0,0 +1,6 @@ +# This is just a reminder. You don't use it. + +# import Pkg; Pkg.add("DocumenterTools") +using DocumenterTools + +DocumenterTools.genkeys(user="ohno", repo="FiniteDifferenceMatrices.jl") diff --git a/dev/docs.jl b/dev/docs.jl new file mode 100644 index 0000000..69f5717 --- /dev/null +++ b/dev/docs.jl @@ -0,0 +1,3 @@ +# Please run `include("./dev/docs.jl")` on RELP. + +run(`julia --project=docs/ -e 'using Pkg; Pkg.activate("./"); cd("docs"); include("make.jl")'`) diff --git a/dev/pkg.jl b/dev/pkg.jl new file mode 100644 index 0000000..61ca829 --- /dev/null +++ b/dev/pkg.jl @@ -0,0 +1,29 @@ +# This is just a reminder. You don't use it. + +# import Pkg; Pkg.add("PkgTemplates") +using PkgTemplates + +t = Template(; + user = "ohno", + authors = ["Shuhei Ohno"], + dir = pwd(), + julia = v"1", + plugins = [ + License(; name = "MIT"), + ProjectFile(; version=v"0.0.1"), + Git(; manifest = false, ssh = true), + GitHubActions(; + extra_versions = ["1.7"] + ), + Documenter{GitHubActions}(), + Readme(; + inline_badges = true, + badge_order = DataType[ + GitHubActions, + Documenter{GitHubActions}, + ], + ), + ], +) + +generate("FiniteDifferenceMatrices", t) diff --git a/dev/revice.jl b/dev/revice.jl new file mode 100644 index 0000000..8cf8422 --- /dev/null +++ b/dev/revice.jl @@ -0,0 +1,10 @@ +# Please run `include("./dev/revice.jl")` on RELP + +# import Pkg; Pkg.add("Revise") +using Revise + +dir = dirname(@__FILE__) * "/../" +cd(dir) +@show pwd() +import Pkg; Pkg.activate(".") +using FiniteDifferenceMatrices diff --git a/dev/test.jl b/dev/test.jl new file mode 100644 index 0000000..81d422d --- /dev/null +++ b/dev/test.jl @@ -0,0 +1,9 @@ +# Please run `include("./dev/test.jl")` on RELP. + +dir = dirname(@__FILE__) * "/../" +cd(dir) +@show pwd() + +using Pkg +Pkg.activate(dir) +Pkg.test() diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 0000000..fe1eb03 --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,3 @@ +[deps] +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +FiniteDifferenceMatrices = "a7a66f33-e7b8-47af-b618-f9b5bea05f3d" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..c91cfe3 --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,24 @@ +using FiniteDifferenceMatrices +using Documenter + +DocMeta.setdocmeta!(FiniteDifferenceMatrices, :DocTestSetup, :(using FiniteDifferenceMatrices); recursive=true) + +makedocs(; + modules=[FiniteDifferenceMatrices], + authors="Shuhei Ohno", + sitename="FiniteDifferenceMatrices.jl", + format=Documenter.HTML(; + canonical="https://ohno.github.io/FiniteDifferenceMatrices.jl", + edit_link="main", + assets=String[], + ), + pages=[ + "Home" => "index.md", + "API reference" => "API.md", + ], +) + +deploydocs(; + repo="github.com/ohno/FiniteDifferenceMatrices.jl", + devbranch="main", +) diff --git a/docs/src/API.md b/docs/src/API.md new file mode 100644 index 0000000..de9e48a --- /dev/null +++ b/docs/src/API.md @@ -0,0 +1,12 @@ +```@meta +CurrentModule = FiniteDifferenceMatrices +``` + +# API reference + +```@index +``` + +```@autodocs +Modules = [FiniteDifferenceMatrices] +``` diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 0000000..997428d --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,72 @@ +```@meta +CurrentModule = FiniteDifferenceMatrices +``` + +# FiniteDifferenceMatrices.jl + +[![Build Status](https://github.com/ohno/FiniteDifferenceMatrices.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/ohno/FiniteDifferenceMatrices.jl/actions/workflows/CI.yml?query=branch%3Amain) +[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://ohno.github.io/FiniteDifferenceMatrices.jl/stable/) +[![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://ohno.github.io/FiniteDifferenceMatrices.jl/dev/) + +A Julia package for discrete approximations of differential operators + +## Install + +Run the following code on the to install this package. +```julia +import Pkg; Pkg.add(url="https://github.com/ohno/FiniteDifferenceMatrices.jl.git") +``` + +## Usage + +Run the following code before each use. +```@example index +using FiniteDifferenceMatrices +``` + +A central finite difference of the second-order derivative is +```math +\frac{\mathrm{d}^{2}f(x)}{\mathrm{d} x^{2}} = \frac{f(x+\Delta x) - 2f(x) + f(x-\Delta x)}{\Delta x^{2}} + O(\Delta x^{2}). +``` +```@repl index +fdcoefficient(n=2, m=2, d=:c) +``` + +A discrete approximation of the second-order differential operator is +```math +\frac{\mathrm{d}^2}{\mathrm{d} x^2} +\simeq +\frac{1}{\Delta x^2} +\left(\begin{array}{ccccccc} + -2 & 1 & 0 & \ldots & 0 & 0 & 0 \\ + 1 & -2 & 1 & \ldots & 0 & 0 & 0 \\ + 0 & 1 & -2 & \ldots & 0 & 0 & 0 \\ + \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ + 0 & 0 & 0 & \ldots & -2 & 1 & 0 \\ + 0 & 0 & 0 & \ldots & 1 & -2 & 1 \\ + 0 & 0 & 0 & \ldots & 0 & 1 & -2 +\end{array}\right). +``` +```@repl index +fdmatrix(5, n=2, m=2, d=:c, h=1//1) +``` + +Please see the [API reference](./API.md) for details and more examples. + +## Developer's Guide + +There are several tools for developers. + +```sh +git clone https://github.com/ohno/FiniteDifferenceMatrices.jl.git +cd FiniteDifferenceMatrices.jl +julia +julia> include("dev/revice.jl") +julia> include("dev/test.jl") +julia> include("dev/docs.jl") +``` + +## API reference + +```@index +``` diff --git a/src/FiniteDifferenceMatrices.jl b/src/FiniteDifferenceMatrices.jl new file mode 100644 index 0000000..46156b8 --- /dev/null +++ b/src/FiniteDifferenceMatrices.jl @@ -0,0 +1,302 @@ +module FiniteDifferenceMatrices + +# export + +export fdcoefficient, fdvalue, fdmatrix + +# package + +using SparseArrays + +# implementation + +@doc raw""" +`Fornberg1988[d,n,m]` + +This is a foolish implementation of the tables of [B. Fornberg, _Math. Comp._ **51** 699-706 (1988)](https://doi.org/10.1090/S0025-5718-1988-0935077-0) and [Wikipedia](https://en.wikipedia.org/wiki/Finite_difference_coefficient) by hand. This was coded to test the `fdcoefficient()`. + +## Examples +```jldoctest +julia> FiniteDifferenceMatrices.Fornberg1988[:c,1,2] +1×3 Matrix{Rational{Int64}}: + -1//2 0//1 1//2 + +julia> FiniteDifferenceMatrices.Fornberg1988[:c,2,2] +1×3 Matrix{Rational{Int64}}: + 1//1 -2//1 1//1 +``` +""" +Fornberg1988 = Dict( + # Central + (:c, 1, 2) => [ -1//2 0//1 1//2 ], + (:c, 1, 4) => [ 1//12 -2//3 0//1 2//3 -1//12 ], + (:c, 1, 6) => [ -1//60 3//20 -3//4 0//1 3//4 -3//20 1//60 ], + (:c, 1, 8) => [ 1//280 -4//105 1//5 -4//5 0//1 4//5 -1//5 4//105 -1//280 ], + (:c, 2, 2) => [ 1//1 -2//1 1//1 ], + (:c, 2, 4) => [ -1//12 4//3 -5//2 4//3 -1//12 ], + (:c, 2, 6) => [ 1//90 -3//20 3//2 -49//18 3//2 -3//20 1//90 ], + (:c, 2, 8) => [ -1//560 8//315 -1//5 8//5 -205//72 8//5 -1//5 8//315 -1//560 ], + (:c, 3, 2) => [ -1//2 1//1 0//1 -1//1 1//2 ], + (:c, 3, 4) => [ 1//8 -1//1 13//8 0//1 -13//8 1//1 -1//8 ], + (:c, 3, 6) => [ -7//240 3//10 -169//120 61//30 0//1 -61//30 169//120 -3//10 7//240 ], + (:c, 4, 2) => [ 1//1 -4//1 6//1 -4//1 1//1 ], + (:c, 4, 4) => [ -1//6 2//1 -13//2 28//3 -13//2 2//1 -1//6 ], + (:c, 4, 6) => [ 7//240 -2//5 169//60 -122//15 91//8 -122//15 169//60 -2//5 7//240 ], + (:c, 5, 2) => [ -1//2 2//1 -5//2 0//1 5//2 -2//1 1//2 ], + (:c, 5, 4) => [ 1//6 -3//2 13//3 -29//6 0//1 29//6 -13//3 3//2 -1//6 ], + (:c, 5, 6) => [-13//288 19//36 -87//32 13//2 -323//48 0//1 323//48 -13//2 87//32 -19//36 13//288], + (:c, 6, 2) => [ 1//1 -6//1 15//1 -20//1 15//1 -6//1 1//1 ], + (:c, 6, 4) => [ -1//4 3//1 -13//1 29//1 -75//2 29//1 -13//1 3//1 -1//4 ], + (:c, 6, 6) => [ 13//240 -19//24 87//16 -39//2 323//8 -1023//20 323//8 -39//2 87//16 -19//24 13//240], + # Forward + (:f, 1, 1) => [ -1//1 1//1 ], + (:f, 1, 2) => [ -3//2 2//1 -1//2 ], + (:f, 1, 3) => [ -11//6 3//1 -3//2 1//3 ], + (:f, 1, 4) => [ -25//12 4//1 -3//1 4//3 -1//4 ], + (:f, 1, 5) => [ -137//60 5//1 -5//1 10//3 -5//4 1//5 ], + (:f, 1, 6) => [ -49//20 6//1 -15//2 20//3 -15//4 6//5 -1//6 ], + (:f, 1, 7) => [ -363//140 7//1 -21//2 35//3 -35//4 21//5 -7//6 1//7 ], + (:f, 1, 8) => [ -761//280 8//1 -14//1 56//3 -35//2 56//5 -14//3 8//7 -1//8 ], + (:f, 2, 1) => [ 1//1 -2//1 1//1 ], + (:f, 2, 2) => [ 2//1 -5//1 4//1 -1//1 ], + (:f, 2, 3) => [ 35//12 -26//3 19//2 -14//3 11//12 ], + (:f, 2, 4) => [ 15//4 -77//6 107//6 -13//1 61//12 -5//6 ], + (:f, 2, 5) => [ 203//45 -87//5 117//4 -254//9 33//2 -27//5 137//180 ], + (:f, 2, 6) => [ 469//90 -223//10 879//20 -949//18 41//1 -201//10 1019//180 -7//10 ], + (:f, 2, 7) => [29531//5040 -962//35 621//10 -4006//45 691//8 -282//5 2143//90 -206//35 363//560], + (:f, 3, 1) => [ -1//1 3//1 -3//1 1//1 ], + (:f, 3, 2) => [ -5//2 9//1 -12//1 7//1 -3//2 ], + (:f, 3, 3) => [ -17//4 71//4 -59//2 49//2 -41//4 7//4 ], + (:f, 3, 4) => [ -49//8 29//1 -461//8 62//1 -307//8 13//1 -15//8 ], + (:f, 3, 5) => [ -967//120 638//15 -3929//40 389//3 -2545//24 268//5 -1849//120 29//15 ], + (:f, 3, 6) => [ -801//80 349//6 -18353//120 2391//10 -1457//6 4891//30 -561//8 527//30 -469//240], + (:f, 4, 1) => [ 1//1 -4//1 6//1 -4//1 1//1 ], + (:f, 4, 2) => [ 3//1 -14//1 26//1 -24//1 11//1 -2//1 ], + (:f, 4, 3) => [ 35//6 -31//1 137//2 -242//3 107//2 -19//1 17//6 ], + (:f, 4, 4) => [ 28//3 -111//2 142//1 -1219//6 176//1 -185//2 82//3 -7//2 ], + (:f, 4, 5) => [ 1069//80 -1316//15 15289//60 -2144//5 10993//24 -4772//15 2803//20 -536//15 967//240], + # Backward + (:b, 1, 1) => [-1//1 1//1 ], + (:b, 1, 2) => [ 1//2 -2//1 3//2 ], + (:b, 1, 3) => [-1//3 3//2 -3//1 11//6 ], + (:b, 2, 1) => [ 1//1 -2//1 1//1 ], + (:b, 2, 2) => [-1//1 4//1 -5//1 2//1 ], + (:b, 3, 1) => [-1//1 3//1 -3//1 1//1 ], + (:b, 3, 2) => [ 3//2 -7//1 12//1 -9//1 5//2 ], + (:b, 4, 1) => [ 1//1 -4//1 6//1 -4//1 1//1 ], + (:b, 4, 2) => [-2//1 11//1 -24//1 26//1 -14//1 3//1], +) + +@doc raw""" +`fdcoefficient(; n::Int=1, m::Int=2, d=:c)` + +This function returns a `Dict` of the finite difference coefficients ``c_i`` of + +```math +\frac{\mathrm{d}^n f(x)}{\mathrm{d}x^n} = \frac{1}{\Delta x^n} \sum_{i} c_i f(x+i\Delta x) + O(\Delta x^m). +``` + +This implementation is based on [a post on discourse](https://discourse.julialang.org/t/generating-finite-difference-stencils/85876/5) by [@stevengj](https://discourse.julialang.org/u/stevengj/summary) and this function is tested to return results equivalent to [B. Fornberg, _Math. Comp._ **51** 699-706 (1988)](https://doi.org/10.1090/S0025-5718-1988-0935077-0). + +| Arguments | Description | +| --- | --- | +| `n` | order of derivative ``n`` | +| `m` | order of accuracy ``m`` | +| `d` | direction, `:c` central, `:f` forward, `:b` backward | + +## Examples + +The coefficients of the central, ``n=1`` and ``m=2`` differences are ``c_{-1} = -1/2, c_{0} = 0, c_{1} = 1/2``. +```math +\frac{\mathrm{d}f(x)}{\mathrm{d} x} = \frac{f(x+\Delta x) - f(x-\Delta x)}{2\Delta x} + O(\Delta x^{2}) +``` + +```jldoctest +julia> fdcoefficient(n=1, m=2, d=:c) +Dict{Int64, Rational{Int64}} with 3 entries: + 0 => 0//1 + -1 => -1//2 + 1 => 1//2 +``` + +The coefficients of the central, ``n=1`` and ``m=1`` differences are ``c_{0} = -1, c_{1} = 1``. + +```math +\frac{\mathrm{d}f(x)}{\mathrm{d} x} = \frac{f(x+\Delta x) - f(x)}{\Delta x} + O(\Delta x) +``` + +```jldoctest +julia> fdcoefficient(n=1, m=1, d=:f) +Dict{Int64, Rational{Int64}} with 2 entries: + 0 => -1//1 + 1 => 1//1 +``` + +The coefficients of the central, ``n=2`` and ``m=2`` differences are ``c_{-1} = 1, c_{0} = -2, c_{1} = 1``. +```math +\frac{\mathrm{d}^{2}f(x)}{\mathrm{d} x^{2}} = \frac{f(x+\Delta x) - 2f(x) + f(x-\Delta x)}{\Delta x^{2}} + O(\Delta x^{2}) +``` + +```jldoctest +julia> fdcoefficient(n=2, m=2, d=:c) +Dict{Int64, Rational{Int64}} with 3 entries: + 0 => -2//1 + -1 => 1//1 + 1 => 1//1 +``` +""" +function fdcoefficient(; n::Int=1, m::Int=2, d=:c) + # `:central` central, `dⁿf/dxⁿ = [f(x+lh) + ... + f(x-lh)] / hⁿ + O(hᵐ), l = Int(m/2) + Int(ceil(n/2)) - 1` + # `:forward` forward, `dⁿf/dxⁿ = [f(x) + ... + f(x+mh)] / hⁿ + O(hᵐ)` + # `:backward` backward, `dⁿf/dxⁿ = [f(x-mh) + ... + f(x) ] / hⁿ + O(hᵐ)` + n ≥ 1 || throw(ArgumentError("n = $n must satisfy n ≥ 1.")) + m ≥ 1 || throw(ArgumentError("m = $m must satisfy m ≥ 1.")) + if d == :c + iseven(m) || throw(ArgumentError("m = $m must be even number for central finite difference.")) + end + I = NaN + I = d == :c ? (-(Int(m/2)+Int(ceil(n/2))-1):(Int(m/2)+Int(ceil(n/2))-1)) : I + I = d == :f ? (0:n+m-1) : I + I = d == :b ? (-n-m+1:0) : I + I₀ = 0//1 + ℓ = 0:length(I)-1 + n in ℓ || throw(ArgumentError("n $n ∉ $ℓ")) + A = @. (I' - I₀)^ℓ / factorial(ℓ) + C = A \ (ℓ .== n) + return Dict(I[j] => C[j] for j in keys(I)) +end + +@doc raw""" +`fdvalue(f, a; n::Int=1, m::Int=2, d=:c, h=0.1)` + +This function calculates the differential coefficient $f^{(n)}(a)$, a value of the derivative function $f^{(n)}=\frac{\mathrm{d}^n f}{\mathrm{d}x^n}$ at the point $a$. + +```math +\frac{\mathrm{d}^n f}{\mathrm{d}x^n} (a) = \frac{1}{h^n} \sum_{i} c_i f(a+ih) + O(\Delta x^m). +``` + +## Examples + +```jldoctest +julia> fdvalue(x -> x^2, 1.0) +2.0000000000000004 + +julia> (x -> 2*x)(1.0) +2.0 + +julia> fdvalue(sin, 0.0) +0.9983341664682815 + +julia> fdvalue(sin, 0.0, m=4) +0.9999966706326067 + +julia> fdvalue(sin, 0.0, m=6) +0.9999999928710186 + +julia> cos(0.0) +1.0 +``` +""" +function fdvalue(f, a; n::Int=1, m::Int=2, d=:c, h=0.1) + C = fdcoefficient(n=n, m=m, d=d) + sum(C[i]*f(a+i*h) for i in keys(C)) / h^n +end + +@doc raw""" +`fdmatrix(N::Int; n::Int=1, m::Int=2, d=:c, h=0.1)` + +This function returns a discrete approximation of a differential operator as a [sparse matrix](https://github.com/JuliaSparse/SparseArrays.jl). + +| Arguments | Description | +| --- | --- | +| `n` | order of derivative ``n`` | +| `m` | order of accuracy ``m`` | +| `d` | direction, `:c` central, `:f` forward, `:b` backward | +| `h` | grid spacing ``\Delta x`` | + +## Examples + +The central, ``n=1`` and ``m=2`` discrete approximation of the differential operator is +```math +\frac{\mathrm{d}}{\mathrm{d} x} +\simeq +\frac{1}{2\Delta x} +\left(\begin{array}{ccccccc} + 0 & 1 & 0 & \ldots & 0 & 0 & 0 \\ + -1 & 0 & 1 & \ldots & 0 & 0 & 0 \\ + 0 & -1 & 0 & \ldots & 0 & 0 & 0 \\ + \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ + 0 & 0 & 0 & \ldots & 0 & 1 & 0 \\ + 0 & 0 & 0 & \ldots & -1 & 0 & 1 \\ + 0 & 0 & 0 & \ldots & 0 & -1 & 0 +\end{array}\right). +``` + +```jldoctest +julia> fdmatrix(5, n=1, m=2, d=:c, h=1//1) +5×5 SparseArrays.SparseMatrixCSC{Rational{Int64}, Int64} with 8 stored entries: + ⋅ 1//2 ⋅ ⋅ ⋅ + -1//2 ⋅ 1//2 ⋅ ⋅ + ⋅ -1//2 ⋅ 1//2 ⋅ + ⋅ ⋅ -1//2 ⋅ 1//2 + ⋅ ⋅ ⋅ -1//2 ⋅ +``` + +The central, ``n=1`` and ``m=1`` discrete approximation of the differential operator is +```math +\frac{\mathrm{d}}{\mathrm{d} x} +\simeq +\frac{1}{\Delta x} +\left(\begin{array}{ccccccc} + -1 & 1 & 0 & \ldots & 0 & 0 & 0 \\ + 0 & -1 & 1 & \ldots & 0 & 0 & 0 \\ + 0 & 0 & -1 & \ldots & 0 & 0 & 0 \\ + \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ + 0 & 0 & 0 & \ldots & -1 & 1 & 0 \\ + 0 & 0 & 0 & \ldots & 0 & -1 & 1 \\ + 0 & 0 & 0 & \ldots & 0 & 0 & -1 +\end{array}\right). +``` + +```jldoctest +julia> fdmatrix(5, n=1, m=1, d=:f, h=1//1) +5×5 SparseArrays.SparseMatrixCSC{Rational{Int64}, Int64} with 9 stored entries: + -1//1 1//1 ⋅ ⋅ ⋅ + ⋅ -1//1 1//1 ⋅ ⋅ + ⋅ ⋅ -1//1 1//1 ⋅ + ⋅ ⋅ ⋅ -1//1 1//1 + ⋅ ⋅ ⋅ ⋅ -1//1 +``` + +The central, ``n=2`` and ``m=2`` discrete approximation of the differential operator is +```math +\frac{\mathrm{d}^2}{\mathrm{d} x^2} +\simeq +\frac{1}{\Delta x^2} +\left(\begin{array}{ccccccc} + -2 & 1 & 0 & \ldots & 0 & 0 & 0 \\ + 1 & -2 & 1 & \ldots & 0 & 0 & 0 \\ + 0 & 1 & -2 & \ldots & 0 & 0 & 0 \\ + \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ + 0 & 0 & 0 & \ldots & -2 & 1 & 0 \\ + 0 & 0 & 0 & \ldots & 1 & -2 & 1 \\ + 0 & 0 & 0 & \ldots & 0 & 1 & -2 +\end{array}\right). +``` + +```jldoctest +julia> fdmatrix(5, n=2, m=2, d=:c, h=1//1) +5×5 SparseArrays.SparseMatrixCSC{Rational{Int64}, Int64} with 13 stored entries: + -2//1 1//1 ⋅ ⋅ ⋅ + 1//1 -2//1 1//1 ⋅ ⋅ + ⋅ 1//1 -2//1 1//1 ⋅ + ⋅ ⋅ 1//1 -2//1 1//1 + ⋅ ⋅ ⋅ 1//1 -2//1 +``` +""" +function fdmatrix(N::Int; n::Int=1, m::Int=2, d=:c, h=0.1) + C = fdcoefficient(n=n, m=m, d=d) + return h^(-n) * spdiagm(Dict(i => fill(C[i],N-abs(i)) for i in keys(C))...) +end + +end diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..8e3c52b --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1,27 @@ +using FiniteDifferenceMatrices +using Test + +@testset "FiniteDifferenceMatrices.jl" begin + + for k in keys(FiniteDifferenceMatrices.Fornberg1988) + d = k[1] + n = k[2] + m = k[3] + I = NaN + I = d == :c ? (-(Int(m/2)+Int(ceil(n/2))-1):(Int(m/2)+Int(ceil(n/2))-1)) : I + I = d == :f ? (0:n+m-1) : I + I = d == :b ? (-n-m+1:0) : I + C = fdcoefficient(n=n, m=m, d=d) + println("d = :$d, n = $n, m = $m") + for j in 1:length(I) + i = I[j] + Fornberg = FiniteDifferenceMatrices.Fornberg1988[k][j] + ThisWork = C[i] + acceptance = isequal(Fornberg, ThisWork) + println(acceptance ? "✔" : "✗", " i = $i $Fornberg == $ThisWork") + @test acceptance + end + println() + end + +end