Skip to content

Commit 2ddba42

Browse files
Anonymousoverleaf
Anonymous
authored and
overleaf
committed
Project created
0 parents  commit 2ddba42

10 files changed

+1143
-0
lines changed

‎.gitignore

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
*.md.pdf
2+
*-cache.yaml
3+
*.aux
4+
*.fls
5+
*.log
6+
*.pdf
7+
*.gz
8+
_minted*
9+
*.fdb_latexmk
10+
*.bbl
11+
*.blg

‎BellmanDiffusion.md

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
# Solving the Bellman Equation with a Simple Univariate Diffusion
2+
## Setup
3+
Take the stochastic process
4+
$$
5+
d x_t = \mu(x_t)dt + \sigma(x_t)d W_t
6+
$$
7+
where $W_t$ is Brownian motion and reflecting barriers at $x \in (\underline{x}, \bar{x})$ where $-\infty < \underline{x} < \bar{x} < \infty$.
8+
9+
The partial differential operator (infinitesimal generator) associated with the stochastic process is
10+
$$
11+
\mathcal{A} \equiv \mu(x)\partial_x + \frac{\sigma(x)^2}{2}\partial_{xx}
12+
$$
13+
14+
Then, if the payoff in state $x$ is $u(x)$, and payoffs are discounted at rate $\rho$, then the Bellman equation is,
15+
$$
16+
\rho v(x) = u(x) + \mathcal{A}v(x)
17+
$$
18+
With boundary values $v'(\underline{x}) = 0$ and $v'(\bar{x}) = 0$
19+
20+
## Solving the discretized problem
21+
Create a grid on $x$ where $i = 1, \ldots I$ and and $x_1 = \underline{x} = x_1, x2, \ldots, x_I$, define $v \in \mathbb{R}^I$ such that $v_i = v(x_i)$, and finally $u \in \mathbb{R}^I$ such that $u_i = u(x_i)$.
22+
23+
Let the upwind discretized $\mathcal{A}$, subject to the boundary conditions, be $A$, then the solution to the bellman equation is the solution to the following linear system of equations,
24+
$$
25+
\rho v = u + A v
26+
$$
27+
28+
```julia
29+
f(x) = x
30+
```

‎DiffEqOperators.md

Lines changed: 138 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,138 @@
1+
# Basics of DiffEqOperators
2+
3+
In this tutorial we will explore the basic functionalities of PDEOperator which is used to obtain the discretizations of PDEs of appropriate derivative and approximation order.
4+
5+
So an operator API is as follows:-
6+
7+
A = DerivativeOperator{T}
8+
(
9+
derivative_order,
10+
approximation_order,
11+
grid_step_size,
12+
grid_size,
13+
:LBC,
14+
:RBC;
15+
BC=(LBV, RBV)
16+
);
17+
Currently we support the `Dirichlet 0/1`, `Neumann 0/1`, `periodic` and `Robin` boundary conditions.
18+
19+
Taking a specific example
20+
21+
A = DerivativeOperator{Float64}(2,2,1/99,10,:Dirichlet,:Dirichlet; BC=(u[1],u[end]))
22+
23+
this is the time independent Dirichlet BC. You can also specify a time dependent Dirichlet BC as follows:-
24+
25+
A = DerivativeOperator{Float64}(2,2,1/99,10,:Dirichlet,:Dirichlet; bndry_fn=(t->(u[1]*cos(t)),u[end]))
26+
27+
We have generated an operator which produces the 2nd order approximation of the Laplacian. We can checkout the stencil as follows:-
28+
29+
julia> A.stencil_coefs
30+
3-element SVector{3,Float64}:
31+
1.0
32+
-2.0
33+
1.0
34+
35+
We can get the linear operator as a matrix as follows:-
36+
37+
julia> full(A)
38+
10×10 Array{Float64,2}:
39+
-2.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
40+
1.0 -2.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
41+
0.0 1.0 -2.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
42+
0.0 0.0 1.0 -2.0 1.0 0.0 0.0 0.0 0.0 0.0
43+
0.0 0.0 0.0 1.0 -2.0 1.0 0.0 0.0 0.0 0.0
44+
0.0 0.0 0.0 0.0 1.0 -2.0 1.0 0.0 0.0 0.0
45+
0.0 0.0 0.0 0.0 0.0 1.0 -2.0 1.0 0.0 0.0
46+
0.0 0.0 0.0 0.0 0.0 0.0 1.0 -2.0 1.0 0.0
47+
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 -2.0 1.0
48+
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 -2.0
49+
50+
Note that we **don't** need to define the `BC` only for `:D0` and `:periodic` boundary conditions so you can ignore it.
51+
52+
53+
Now coming to the main functionality of DiffEqOperators ie. taking finite difference discretizations of functions.
54+
55+
julia> x = collect(0 : 1/99 : 1);
56+
julia> u0 = x.^2 -x;
57+
julia> res = A*u0
58+
100-element Array{Float64,1}:
59+
-98.0
60+
2.0
61+
2.0
62+
2.0
63+
2.0
64+
2.0
65+
2.0
66+
⋮
67+
2.0
68+
2.0
69+
2.0
70+
2.0
71+
2.0
72+
2.0
73+
-98.0
74+
75+
The derivative values at the boundaries are in accordance with the `Dirichlet` boundary condition.
76+
77+
You can also take derivatives of matrices using `A*M` or `M*A` where the order of multiplication decides the axis along which we want to take derivatives.
78+
79+
julia> xarr = linspace(0,1,51)
80+
julia> yarr = linspace(0,1,101)
81+
julia> dx = xarr[2]-xarr[1]
82+
julia> dy = yarr[2]-yarr[1]
83+
julia> F = [x^2+y for x = xarr, y = yarr]
84+
julia> A = DerivativeOperator{Float64}(2,2,dx,length(yarr),:None,:None)
85+
julia> B = DerivativeOperator{Float64}(2,2,dy,length(yarr),:None,:None)
86+
87+
julia> # A*F calculates derivatives along the x axis ie. keeping y constant
88+
julia> A*F
89+
51×101 Array{Float64,2}:
90+
2.0 2.0 2.0 2.0 2.0 2.0 2.0 … 2.0 2.0 2.0 2.0 2.0 2.0 2.0
91+
2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
92+
2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
93+
2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
94+
2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
95+
2.0 2.0 2.0 2.0 2.0 2.0 2.0 … 2.0 2.0 2.0 2.0 2.0 2.0 2.0
96+
2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
97+
2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
98+
2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
99+
⋮ ⋮ ⋱ ⋮ ⋮
100+
2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
101+
2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
102+
2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
103+
2.0 2.0 2.0 2.0 2.0 2.0 2.0 … 2.0 2.0 2.0 2.0 2.0 2.0 2.0
104+
2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
105+
2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
106+
2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
107+
2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0
108+
2.0 2.0 2.0 2.0 2.0 2.0 2.0 … 2.0 2.0 2.0 2.0 2.0 2.0 2.0
109+
110+
julia> # F*B calculates derivatives along the y axis ie. keeping x constant
111+
julia> F*B
112+
101×51 Array{Float64,2}:
113+
0.0 1.04083e-13 6.93889e-14 … 2.22045e-12 2.22045e-12
114+
0.0 3.46945e-14 0.0 0.0 0.0
115+
0.0 -3.46945e-14 -6.93889e-14 0.0 0.0
116+
6.93889e-14 0.0 6.93889e-14 0.0 0.0
117+
0.0 0.0 0.0 0.0 0.0
118+
-6.93889e-14 -6.93889e-14 -6.93889e-14 … 0.0 0.0
119+
1.38778e-13 1.38778e-13 1.38778e-13 0.0 0.0
120+
-1.38778e-13 -1.38778e-13 -2.77556e-13 0.0 0.0
121+
0.0 0.0 0.0 0.0 0.0
122+
⋮ ⋱ ⋮
123+
0.0 0.0 0.0 4.44089e-12 4.44089e-12
124+
-1.11022e-12 -1.11022e-12 -1.11022e-12 -4.44089e-12 -4.44089e-12
125+
1.11022e-12 1.11022e-12 1.11022e-12 2.22045e-12 2.22045e-12
126+
0.0 0.0 0.0 … 0.0 0.0
127+
0.0 0.0 0.0 0.0 0.0
128+
0.0 0.0 0.0 0.0 0.0
129+
0.0 0.0 0.0 0.0 0.0
130+
0.0 0.0 0.0 0.0 0.0
131+
0.0 0.0 0.0 … 8.88178e-12 8.88178e-12
132+
133+
134+
135+
**Note:** Please take care that the boundary values passed to the operator match the initial boundary conditions. The operator with the boundary condition is meant to enforce the boundary condition rather bring the boundaries to that state. ~~Right now we support only **constant** boundaries conditions, time dependent conditions will supported in later versions.~~
136+
Support for time dependent Dirichlet BC has been added.
137+
138+
**Note:** If you want to parallelize the operation of PDEOperator, please start Julia by specifying the number of threads using `export JULIA_NUM_THREADS=<desired number of threads>`

‎HeatEquation.md

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,21 @@
1+
# Solving the Heat Equation using DiffEqOperators
2+
3+
In this tutorial we will solve the famous heat equation using the explicit discretization on a 2D `space x time` grid. The heat equation is:-
4+
5+
$$\frac{\partial u}{\partial t} - \frac{{\partial}^2 u}{\partial x^2} = 0$$
6+
7+
For this example we consider a Dirichlet boundary condition with the initial distribution being parabolic. Since we have fixed the value at boundaries (in this case equal), after a long time we expect the 1D rod to be heated in a linear manner.
8+
9+
julia> using DiffEqOperators, DifferentialEquations, Plots
10+
julia> x = collect(-pi : 2pi/511 : pi);
11+
julia> u0 = -(x - 0.5).^2 + 1/12;
12+
julia> A = DerivativeOperator{Float64}(2,2,2pi/511,512,:Dirichlet,:Dirichlet;BC=(u0[1],u0[end]));
13+
14+
Now solving equation as an ODE we have:-
15+
16+
julia> prob1 = ODEProblem(A, u0, (0.,10.));
17+
julia> sol1 = solve(prob1, dense=false, tstops=0:0.01:10);
18+
# try to plot the solution at different time points using
19+
julia> plot(x, [sol1(i) for i in 0:1:10])
20+
21+
**Note:** Many a times the solver may inform you that the solution is unstable. This problem is usually handled by changing the solver algorithm. Your best bet might be the `CVODE_BDE()` algorithm from the OrdinaryDiffEq.jl suite.

‎KdV.md

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
# Solving the Heat Equation using DiffEqOperators
2+
3+
In this tutorial we will try to solve the famous **KdV equation** which describes the motion of waves on shallow water surfaces.
4+
The equation is commonly written as
5+
6+
$$\frac{\partial u}{\partial t} + \frac{{\partial}^3 u}{\partial t^3} - 6*u*\frac{\partial u}{\partial t} = 0$.
7+
8+
Lets consider the cosine wave as the initial waveform and evolve it using the equation with a [periodic boundary condition](https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.15.240). This example is taken from [here](https://en.wikipedia.org/wiki/Korteweg%E2%80%93de_Vries_equation).
9+
10+
using DiffEqOperators, DifferentialEquations, Plots
11+
x = collect(0 : 1/99 : 2);
12+
u0 = cos.(π*x);
13+
du3 = zeros(u0); # a container array
14+
function KdV(t, u, du)
15+
C(t,u,du3)
16+
A(t, u, du)
17+
copy!(du, -u.*du .- 0.022^2.*du3)
18+
end
19+
20+
Now defining our DiffEqOperators
21+
```
22+
A = DerivativeOperator{Float64}(1,2,1/99,199,:periodic,:periodic);
23+
C = DerivativeOperator{Float64}(3,2,1/99,199,:periodic,:periodic);
24+
```
25+
26+
Now call the ODE solver as follows:-
27+
28+
prob1 = ODEProblem(KdV, u0, (0.,5.));
29+
sol1 = solve(prob1, dense=false, tstops=0:0.01:10);
30+
31+
and plot the solutions to see the waveform evolve with time.
32+
**Note:** The waveform being solved for here is non-directional unlike the many waves you might see like traveling solitons. In that case you might need to use the Upwind operators.

‎Telegrapher_Equations.ipynb

Lines changed: 244 additions & 0 deletions
Large diffs are not rendered by default.

0 commit comments

Comments
 (0)