ElementaryPDESolutions

Compute polynomial solutions to certain PDEs given a polynomial source term.

Overview

This package provides functionality for solving

\[ \mathcal{L}P = Q,\]

where Q is a source term of Polynomial type (or a NTuple of polynomials for vector-valued problems), P is the sought polynomial solution, and $\mathcal{L}$ is a certain (linear) constant coefficient differential operator.

A typical use case consists of creating a polynomial Q

julia> using ElementaryPDESolutions;
julia> Q = Polynomial([(1,2)=>2, (3,0)=>-1])2xy² - x³

and calling the desired method to obtain P:

julia> P = solve_helmholtz(Q;k=1)2.0x + 2.0xy² - x³

Note that P can be used as function:

julia> P((0.1,0.2)) # functor interface0.20700000000000002

The following PDEs are currently implemented (see their specific Docstrings for further details):

• solve_laplace
• solve_anisotropic_laplace
• solve_anisotropic_advect
• solve_anisotropic_advect_diffuse
• solve_helmholtz
• solve_bilaplace
• solve_stokes
• solve_brinkman
• solve_elastostatic
• solve_elastodynamics
• solve_maxwell

Coefficient type, precision, and conversions

The type of coefficients in the polynomial solution P is inferred from both the type of the coefficients of the input polynomial Q, and from the parameter types (e.g. the type of k in solve_helmholtz). In the presence of floating point errors, this means that the computed coefficients may be inexact:

julia> using ElementaryPDESolutions;
julia> Q = Polynomial((2,2)=>1)x²y²
julia> P = solve_helmholtz(Q;k=3)0.010973936899862825 - 0.024691358024691357y² - 0.024691358024691357x² + 0.1111111111111111x²y²

The recursive algorithm in solve_helmholtz performs repeated divisions by k²; thus, passing a rational types yields an exact result in this case:

julia> Q = Polynomial((3,2)=>big(1//1))x³y²
julia> P = solve_helmholtz(Q;k=3//1)8//243x - 2//27xy² - 2//81x³ + 1//9x³y²

You can still convert the coefficients back to e.g. a floating point type:

julia> convert(Polynomial{2,Float64},P)0.03292181069958848x - 0.07407407407407407xy² - 0.024691358024691357x³ + 0.1111111111111111x³y²

Alternatively, you can use intervals to obtain rigorous bounds on the coefficients:

julia> using IntervalArithmetic
julia> Q = Polynomial((3,2)=>1)x³y²
julia> P = solve_helmholtz(Q;k=3..3)ERROR: UndefVarError: `..` not defined