# Linear equations

Linear equations such as ordinary and partial differential equations, fractional differential equations and integral equations can be solved using ApproxFun. This is accomplished using `A\b`

where `A`

is an `Operator`

and `b`

is a `Fun`

. As a simple example, consider the equation

\[\mathop{u}'(θ) + c \mathop{u}(θ) = \cos{θ},\]

where we want a solution that is periodic on $[0,2π)$. This can be solved succinctly as follows:

`julia> b = Fun(cos,Fourier());`

`julia> c = 0.1; u = (𝒟+c*I) \ b;`

`julia> u(0.6)`

`0.6407683513722799`

`julia> (c*cos(0.6)+sin(0.6)) / (1+c^2) # exact solution`

`0.6407683513722804`

Recall that `𝒟`

is an alias to `Derivative() == Derivative(UnsetSpace(),1)`

.

As another example, consider the Fredholm integral equation

\[\mathop{u} + \mathop{e}^x \int_{-1}^1 \cos{x} \mathop{u}(x) \mathop{dx} = \cos{\mathop{e}^x}.\]

We can solve this equation as follows:

`julia> Σ = DefiniteIntegral(Chebyshev()); x = Fun();`

`julia> u = (I+exp(x)*Σ[cos(x)]) \ cos(exp(x));`

`julia> u(0.1)`

`0.21864294855628844`

Note that we used the syntax `op[f::Fun]`

, which is a shorthand for `op*Multiplication(f)`

.

## Boundary conditions

Incorporating boundary conditions into differential equations is important so that the equation is well-posed. This is accomplished via combining operators and functionals (i.e., `1 × ∞`

operators). As a simple example, consider the first order initial value problem

\[\begin{gathered} \mathop{u}' = t \mathop{u}, \\ \mathop{u}(0) = 1. \end{gathered}\]

To pose this in ApproxFun, we want to find a `u`

such that `Evaluation(0)*u == 1`

and `(𝒟 - t)*u == 0`

. This is accomplished via:

`julia> t = Fun(0..1);`

`julia> u = [Evaluation(0); 𝒟-t] \ [1;0];`

`julia> u(0)`

`0.9999999999999997`

`julia> norm(u'-t*u)`

`1.1860123917817563e-17`

Behind the scenes, the `Vector{Operator{T}}`

representing the functionals and operators are combined into a single `InterlaceOperator`

.

A common usage is two-point boundary value problems. Consider the singularly perturbed boundary value problem:

\[\begin{gathered} ϵ \mathop{u}'' - x \mathop{u}' + \mathop{u} = 0, \\ \mathop{u}(-1) = 1, \mathop{u}(1) = 2. \end{gathered}\]

This can be solved in ApproxFun via:

`julia> ϵ = 1/70; x = Fun();`

`julia> u = [Evaluation(-1); Evaluation(1); ϵ*𝒟^2-x*𝒟+I] \ [1,2,0];`

`julia> u(0.1)`

`0.049999999999959965`

Note in this case the space is inferred from the variable coefficient `x`

.

This ODE can also be solved using the `Dirichlet`

operator:

`julia> x = Fun();`

`julia> u = [Dirichlet(); 1/70*𝒟^2-x*𝒟+I] \ [[1,2],0];`

`julia> u(0.1)`

`0.049999999999959965`

## Eigenvalue Problems

In analogy to linear algebra, many differential equations may be posed as eigenvalue problems. That is, for some differential operator $\mathop{L}$, there are a family of functions $\mathop{u}_i(x)$ such that

\[\mathop{L} \mathop{u}_i(x) = λ_i \mathop{u}_i(x),\]

where $λ_i$ is the $i^{th}$ eigenvalue of the $L$ and has a corresponding *eigenfunction* $\mathop{u}_i(x)$. A classic eigenvalue problem is known as the quantum harmonic oscillator where

\[\mathop{L} = -\frac{1}{2}\frac{\mathop{d}^2}{\mathop{dx}^2} + \frac{1}{2} x^2,\]

and one demands that $\mathop{u}(∞) = \mathop{u}(-∞) = 0$. Because we expect the solutions to be exponentially suppressed for large $x$, we can approximate this with Dirichlet boundary conditions at a 'reasonably large' $x$ without much difference.

We can express this in ApproxFun as the following:

```
x = Fun(-8..8)
L = -𝒟^2/2 + x^2/2
S = space(x)
B = Dirichlet(S)
λ, v = ApproxFun.eigs(B, L, 500,tolerance=1E-10)
```

Note that boundary conditions must be specified in the call to `eigs`

. Plotting the first $20$ eigenfunctions offset vertically by their eigenvalue, we see

If the solutions are not relatively constant near the boundary then one should push the boundaries further out.

For problems with different contraints or boundary conditions, `B`

can be any zero functional constraint, e.g., `DefiniteIntegral()`

.

## Systems of equations

Systems of equations can be handled by creating a matrix of operators and functionals. For example, we can solve the system

\[\begin{gathered} \mathop{u}'' - \mathop{u} + 2 \mathop{v} = \mathop{e}^x, \\ \mathop{u} + \mathop{v}' + \mathop{v} = \cos{x}, \\ \mathop{u}(-1) = \mathop{u}'(-1) = \mathop{v}(-1) = 0 \end{gathered}\]

using the following code:

`julia> x = Fun(); B = Evaluation(Chebyshev(),-1);`

`julia> A = [B 0; B*𝒟 0; 0 B; 𝒟^2-I 2I; I 𝒟+I];`

`julia> u,v = A \ [0;0;0;exp(x);cos(x)];`

`julia> u(-1),u'(-1),v(-1)`

`(6.938893903907228e-17, -1.1102230246251565e-16, 2.220446049250313e-16)`

`julia> norm(u''-u+2v-exp(x))`

`1.7467412827384798e-16`

`julia> norm(u+v'+v-cos(x))`

`1.3782415611264478e-16`

In this example, the automatic space detection failed and so we needed to specify explicitly that the domain space for `B`

is `Chebyshev()`

.

## QR Factorization

Behind the scenes, `A\b`

where `A`

is an `Operator`

is implemented via an adaptive QR factorization. That is, it is equivalent to `qr(A)\b`

. (There is a subtlety here in space inferring: `A\b`

can use both `A`

and `b`

to determine the domain space, while `qr(A)`

only sees the operator `A`

.) Note that `qr`

adaptively caches a partial QR Factorization as it is applied to different right-hand sides, so the same operator can be inverted much more efficiently in subsequent problems.

## Partial differential equations

Partial differential operators are also supported. Here's an example of solving the Poisson equation with zero boundary conditions:

```
d = Domain(-1..1)^2
x,y = Fun(d)
f = exp.(-10(x+0.3)^2-20(y-0.2)^2) # use broadcasting as exp(f) not implemented in 2D
A = [Dirichlet(d);Δ] # Δ is an alias for Laplacian()
@time u = A \ [zeros(∂(d));f] # 4s for ~3k coefficients
```

Using a QR Factorization reduces the cost of subsequent calls substantially:

```
QR = qr(A)
@time QR \ [zeros(∂(d));f] # 4s
g = exp.(-10(x+0.2)^2-20(y-0.1)^2)
@time QR \ [zeros(∂(d));g] # 0.09s
```

Many PDEs have weak singularities at the corners, in which case it is beneficial to specify a tolerance to reduce the time:

`\(A, [zeros(∂(d));f]; tolerance=1E-6)`

## Nonlinear equations

There is preliminary support for nonlinear equations, via Newton iteration in function space. Here is a simple two-point boundary value problem:

\[\begin{gathered} ϵ \mathop{u}'' + 6(1-x^2) \mathop{u}' + \mathop{u}^2=1, \\ \mathop{u}(-1) = \mathop{u}(1) = 0. \end{gathered}\]

This can be solved as follows:

`julia> ϵ = 1/70; x = Fun();`

`julia> N = u -> [u(-1); u(1); ϵ*u''+6*(1-x^2)*u'+u^2-1];`

`julia> u0 = x; u = newton(N,u0);`

`julia> u(-1), u(1)`

`(-4.440892098500626e-16, 0.0)`

`julia> norm(ϵ*u''+6*(1-x^2)*u'+u^2-1)`

`1.519404143120108e-15`