Skip to content

Commit e23fb8e

Browse files
jishnubjishnub
authored
Revert ylim setting in recipe (#825)
* Fix domain eltype in Poisson PDE example Also, revert ylim setting in recipe * change Poisson domain to Int
1 parent 812c58b commit e23fb8e

18 files changed

+174
-162
lines changed

‎docs/make.jl‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -16,12 +16,12 @@ end
1616

1717
for (example, included) in [
1818
("ODE.jl", ["ODE_BVP.jl", "ODE_increaseprec.jl"]),
19-
("PDE.jl", ["PDE1.jl"]),
20-
("Sampling.jl", String[]),
19+
("PDE.jl", ["PDE_Poisson.jl", "PDE_Helmholtz.jl"]),
20+
("Sampling.jl", ["Sampling1.jl"]),
2121
("Periodic.jl", ["Periodic1.jl"]),
22-
("Eigenvalue.jl", String[]),
23-
("NonlinearBVP.jl", ["NonlinearBVP1.jl"]),
24-
("system_of_eqn.jl", String[])
22+
("Eigenvalue.jl", ["Eigenvalue_standard.jl", "Eigenvalue_symmetric.jl"]),
23+
("NonlinearBVP.jl", ["NonlinearBVP1.jl", "NonlinearBVP2.jl"]),
24+
("system_of_eqn.jl", ["System1.jl"])
2525
]
2626
filename = joinpath(example_dir, example)
2727
Literate.markdown(filename, output_dir, documenter=true,

‎examples/Eigenvalue.jl‎

Lines changed: 4 additions & 71 deletions
Original file line numberDiff line numberDiff line change
@@ -1,42 +1,15 @@
11
# # Eigenvalue problem
22

3-
# ## Standard eigenvalue problem
4-
5-
# In analogy to linear algebra, many differential equations may be posed as eigenvalue problems.
6-
# That is, for some differential operator ``\mathop{L}``, there are a family of functions
7-
# ``\mathop{u}_i(x)`` such that
8-
# ```math
9-
# \mathop{L} \mathop{u}_i(x) = λ_i \mathop{u}_i(x),
10-
# ```
11-
# where ``λ_i`` is the ``i^{th}`` eigenvalue of the ``L`` and has a corresponding
12-
# *eigenfunction* ``\mathop{u}_i(x)``.
13-
# A classic eigenvalue problem is known as the quantum harmonic oscillator where
14-
# ```math
15-
# \mathop{L} = -\frac{1}{2}\frac{\mathop{d}^2}{\mathop{dx}^2} + \frac{1}{2} x^2,
16-
# ```
17-
# and one demands that ``\mathop{u}(∞) = \mathop{u}(-∞) = 0``.
18-
# Because we expect the solutions to be exponentially suppressed for large ``x``,
19-
# we can approximate this with Dirichlet boundary conditions at a 'reasonably large' ``x``
20-
# without much difference.
21-
22-
# We can express this in ApproxFun as the following:
23-
using ApproxFun
24-
using LinearAlgebra
25-
26-
x = Fun(-8..8)
27-
V = x^2/2
28-
L = -𝒟^2/2 + V
29-
S = space(x)
30-
B = Dirichlet(S)
31-
λ, v = ApproxFun.eigs(B, L, 500,tolerance=1E-10);
3+
include("Eigenvalue_standard.jl")
324

335
# Note that boundary conditions must be specified in the call to `eigs`.
346
# Plotting the first ``20`` eigenfunctions offset vertically by their eigenvalue, we see
357

368
import Plots
9+
using LinearAlgebra: norm
3710
p = Plots.plot(V; legend=false, ylim=(-Inf, λ[22]))
3811
for k=1:20
39-
Plots.plot!(real(v[k]/norm(v[k]) + λ[k]))
12+
Plots.plot!(real(v[k]/norm(v[k]) + λ[k]), )
4013
end
4114
p
4215

@@ -46,47 +19,7 @@ p
4619
# For problems with different contraints or boundary conditions,
4720
# `B` can be any zero functional constraint, e.g., `DefiniteIntegral()`.
4821

49-
# ## Self-adjoint Eigenvalue Problem
50-
# Ref:
51-
# [J. L. Aurentz & R. M. Slevinsky (2019), arXiv:1903.08538](https://arxiv.org/abs/1903.08538)
52-
53-
# We solve the confined anharmonic oscillator
54-
# ```math
55-
# \left[-\frac{d^2}{dx^2} + V(x)\right] u = λu,
56-
# ```
57-
# where ``u(\pm 10) = 0``, ``V(x) = ωx^2 + x^4``, and ``ω = 25``.
58-
59-
# Define parameters
60-
ω = 25.0
61-
d = -10..10;
62-
S = Legendre(d) # Equivalently, Ultraspherical(0.5, d)
63-
NS = NormalizedPolynomialSpace(S) # NormalizedLegendre
64-
D1 = Conversion(S, NS)
65-
D2 = Conversion(NS, S);
66-
67-
# Construct the differential operator
68-
V = Fun(x -> ω * x^2 + x^4, S)
69-
L = -Derivative(S, 2) + V;
70-
71-
# Basis recombination to transform to one that satisfies Dirichlet boundary conditions
72-
B = Dirichlet(S)
73-
QS = QuotientSpace(B)
74-
Q = Conversion(QS, S)
75-
R = D1*Q;
76-
77-
# This inversion is computed approximately, such that
78-
# ``\mathrm{C}^{-1} \mathrm{C} ≈ \mathrm{I}`` up to a certain bandwidth
79-
C = Conversion(domainspace(L), rangespace(L))
80-
P = cache(PartialInverseOperator(C, (0, ApproxFun.bandwidth(L, 1) + ApproxFun.bandwidth(R, 1) + ApproxFun.bandwidth(C, 2))));
81-
82-
A = R'D1*P*L*D2*R
83-
B = R'R;
84-
85-
# We impose a cutoff to obtain approximate finite matrix representations
86-
n = 3000
87-
SA = Symmetric(A[1:n,1:n], :L)
88-
SB = Symmetric(B[1:n,1:n], :L)
89-
λ = eigvals(SA, SB)[1:round(Int, 3n/5)];
22+
include("Eigenvalue_symmetric.jl")
9023

9124
# We plot the eigenvalues
9225
import Plots

‎examples/Eigenvalue_standard.jl‎

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,28 @@
1+
# ## Standard eigenvalue problem
2+
3+
# In analogy to linear algebra, many differential equations may be posed as eigenvalue problems.
4+
# That is, for some differential operator ``\mathop{L}``, there are a family of functions
5+
# ``\mathop{u}_i(x)`` such that
6+
# ```math
7+
# \mathop{L} \mathop{u}_i(x) = λ_i \mathop{u}_i(x),
8+
# ```
9+
# where ``λ_i`` is the ``i^{th}`` eigenvalue of the ``L`` and has a corresponding
10+
# *eigenfunction* ``\mathop{u}_i(x)``.
11+
# A classic eigenvalue problem is known as the quantum harmonic oscillator where
12+
# ```math
13+
# \mathop{L} = -\frac{1}{2}\frac{\mathop{d}^2}{\mathop{dx}^2} + \frac{1}{2} x^2,
14+
# ```
15+
# and one demands that ``\mathop{u}(∞) = \mathop{u}(-∞) = 0``.
16+
# Because we expect the solutions to be exponentially suppressed for large ``x``,
17+
# we can approximate this with Dirichlet boundary conditions at a 'reasonably large' ``x``
18+
# without much difference.
19+
20+
# We can express this in ApproxFun as the following:
21+
using ApproxFun
22+
23+
x = Fun(-8..8)
24+
V = x^2/2
25+
L = -𝒟^2/2 + V
26+
S = space(x)
27+
B = Dirichlet(S)
28+
λ, v = ApproxFun.eigs(B, L, 500,tolerance=1E-10);

‎examples/Eigenvalue_symmetric.jl‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
# ## Self-adjoint Eigenvalue Problem
2+
# Ref:
3+
# [J. L. Aurentz & R. M. Slevinsky (2019), arXiv:1903.08538](https://arxiv.org/abs/1903.08538)
4+
5+
# We solve the confined anharmonic oscillator
6+
# ```math
7+
# \left[-\frac{d^2}{dx^2} + V(x)\right] u = λu,
8+
# ```
9+
# where ``u(\pm 10) = 0``, ``V(x) = ωx^2 + x^4``, and ``ω = 25``.
10+
11+
using ApproxFun
12+
using LinearAlgebra
13+
14+
# Define parameters
15+
ω = 25.0
16+
d = -10..10;
17+
S = Legendre(d) # Equivalently, Ultraspherical(0.5, d)
18+
NS = NormalizedPolynomialSpace(S) # NormalizedLegendre
19+
D1 = Conversion(S, NS)
20+
D2 = Conversion(NS, S);
21+
22+
# Construct the differential operator
23+
V = Fun(x -> ω * x^2 + x^4, S)
24+
L = -Derivative(S, 2) + V;
25+
26+
# Basis recombination to transform to one that satisfies Dirichlet boundary conditions
27+
B = Dirichlet(S)
28+
QS = QuotientSpace(B)
29+
Q = Conversion(QS, S)
30+
R = D1*Q;
31+
32+
# This inversion is computed approximately, such that
33+
# ``\mathrm{C}^{-1} \mathrm{C} ≈ \mathrm{I}`` up to a certain bandwidth
34+
C = Conversion(domainspace(L), rangespace(L))
35+
P = cache(PartialInverseOperator(C, (0, ApproxFun.bandwidth(L, 1) + ApproxFun.bandwidth(R, 1) + ApproxFun.bandwidth(C, 2))));
36+
37+
A = R'D1*P*L*D2*R
38+
B = R'R;
39+
40+
# We impose a cutoff to obtain approximate finite matrix representations
41+
n = 3000
42+
SA = Symmetric(A[1:n,1:n], :L)
43+
SB = Symmetric(B[1:n,1:n], :L)
44+
λ = eigvals(SA, SB)[1:round(Int, 3n/5)];

‎examples/NonlinearBVP.jl‎

Lines changed: 2 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -15,23 +15,9 @@ include("NonlinearBVP1.jl")
1515
import Plots
1616
Plots.plot(u; title = "Solution", xlabel="x", ylabel="u(x)", legend=false)
1717

18-
# ## System of nonlinear differential equations
19-
# One can also solve a system of nonlinear ODEs with potentially nonlinear boundary conditions:
20-
N(u1, u2) = [u1'(0) - 0.5*u1(0)*u2(0);
21-
u2'(0) + 1;
22-
u1(1) - 1;
23-
u2(1) - 1;
24-
u1'' + u1*u2;
25-
u2'' - u1*u2]
26-
27-
function nbvpsolver2()
28-
x = Fun(0..1)
29-
u10 = one(x)
30-
u20 = one(x)
31-
newton(N, [u10,u20])
32-
end
33-
u1,u2 = nbvpsolver2();
18+
include("NonlinearBVP2.jl")
3419

3520
# Plot the solutions
21+
import Plots
3622
Plots.plot(u1, label="u1", xlabel="x")
3723
Plots.plot!(u2, label="u2")

‎examples/NonlinearBVP2.jl‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
# ## System of nonlinear differential equations
2+
# One can also solve a system of nonlinear ODEs with potentially nonlinear boundary conditions:
3+
4+
using ApproxFun
5+
using LinearAlgebra
6+
7+
N(u1, u2) = [u1'(0) - 0.5*u1(0)*u2(0);
8+
u2'(0) + 1;
9+
u1(1) - 1;
10+
u2(1) - 1;
11+
u1'' + u1*u2;
12+
u2'' - u1*u2]
13+
14+
function nbvpsolver2()
15+
x = Fun(0..1)
16+
u10 = one(x)
17+
u20 = one(x)
18+
newton(N, [u10,u20])
19+
end
20+
u1,u2 = nbvpsolver2();

‎examples/ODE.jl‎

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -9,8 +9,6 @@
99
# and ``y(b)=\mathrm{Ai}(b)``, where ``\mathrm{Ai}`` represents the Airy function of the
1010
# first kind.
1111

12-
using ApproxFun
13-
1412
include("ODE_BVP.jl")
1513

1614
# We plot the solution
@@ -24,5 +22,5 @@ Plots.plot(u, xlabel="x", ylabel="u(x)", legend=false)
2422
# The following calculates ``e`` to 300 digits by solving the ODE ``u^\prime = u``:
2523

2624
include("ODE_increaseprec.jl")
27-
25+
import Plots
2826
Plots.plot(u; legend=false, xlabel="x", ylabel="u(x)")

‎examples/ODE_BVP.jl‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
1-
using ApproxFun #src
1+
using ApproxFun
2+
using SpecialFunctions
23
# Construct the domain
34
a, b = -20, 10
45
d = a..b;
@@ -16,7 +17,6 @@ B = Dirichlet(d);
1617

1718
# The right hand side for the boundary condition is obtained by evaluating the Airy function
1819
# We use `airyai` from `SpecialFunctions.jl` for this
19-
using SpecialFunctions
2020
B_vals = [airyai(a), airyai(b)];
2121

2222
# The main step, where we solve the differential equation

‎examples/ODE_increaseprec.jl‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,8 @@
1-
using ApproxFun #src
1+
using ApproxFun
22
u = setprecision(1000) do
33
d = BigFloat(0)..BigFloat(1)
44
D = Derivative(d)
55
[ldirichlet(); D-1] \ [1; 0]
6-
end
6+
end;
77
using Test #src
88
@test u(1) ≈ exp(BigFloat(1)) #src

‎examples/PDE.jl‎

Lines changed: 2 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -1,27 +1,6 @@
11
# # Solving PDEs
22

3-
# ## Poisson equation
4-
5-
# We solve the Poisson equation
6-
# ```math
7-
# Δu(x,y) = 0
8-
# ```
9-
# on the rectangle `-1..1 × -1..1`,
10-
# with zero boundary conditions:
11-
12-
using ApproxFun
13-
using LinearAlgebra
14-
15-
d = (-1..1)^2
16-
x,y = Fun(d)
17-
f = exp.(-10(x+0.3)^2-20(y-0.2)^2) # use broadcasting as exp(f) not implemented in 2D
18-
A = [Dirichlet(d); Laplacian()]
19-
u = A \ [zeros(∂(d)); f];
20-
21-
# Using a QR Factorization reduces the cost of subsequent calls substantially
22-
23-
QR = qr(A)
24-
u = QR \ [zeros(∂(d)); f];
3+
include("PDE_Poisson.jl")
254

265
import Plots
276
Plots.plot(u; st=:surface, legend=false)
@@ -31,16 +10,7 @@ Plots.plot(u; st=:surface, legend=false)
3110

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

34-
# ## Helmholtz Equation
35-
36-
# We solve
37-
# ```math
38-
# (Δ + 100)u(x,y) = 0
39-
# ```
40-
# on the rectangle `-1..1 × -1..1`,
41-
# subject to the condition that ``u=1`` on the boundary of the domain.
42-
43-
include("PDE1.jl")
13+
include("PDE_Helmholtz.jl")
4414

4515
# Plot the solution
4616
import Plots

0 commit comments

Comments
 (0)