# Solving Differential Equations

This notebook is less a project for you to work on, much rather a show case of multiple packages, which allow you to solve partial differential equations in Julia. Based on "classic" example in the respective fields I invite you to start playing by yourself. Resources for the respective packages are given.

## Using ModelingToolkit and OrdinaryDiffEq to solve an ODE system

Consider the Lorenz attractor:

$$
\begin{align}
    \dot{x} &= \sigma y - \sigma x\\
    \dot{y} &= \rho x - xz - y\\
    \dot{z} &= xy - \beta z
\end{align}
$$


```julia
using ModelingToolkit
using OrdinaryDiffEq
using Plots

# Define the system to be solved
@parameters t σ ρ β
@variables x(t) y(t) z(t)
D = Differential(t)

eqs = [D(x) ~ σ*y - σ*x,
       D(y) ~ x*ρ - x*z - y,
       D(z) ~ x*y - β*z]
@named lorenz = ODESystem(eqs)

# Define the inital conditions
u0 = [x => 1.0,
      y => 1.0,
      z => 1.0]

# Define the values of the parameters
p  = [σ => 10.0,
      ρ => 28.0,
      β => 8/3]

# Timespan for propagation
tspan = (0.0, 30.0)

# Solve it
prob  = ODEProblem(lorenz, u0, tspan, p, jac=true)
sol   = OrdinaryDiffEq.solve(prob, Tsit5())

# and plot solution
plot(sol, vars=(x,y))
```

##### More details
- https://github.com/SciML/OrdinaryDiffEq.jl
- https://github.com/SciML/DifferentialEquations.jl
- https://github.com/SciML/ModelingToolkit.jl

## Solving PDEs using the finite-element method

Classic example:

$$
\begin{align}
    - \Delta u &= f &&\text{in $\Omega$} \\
    u &= g &&\text{on $\partial \Omega$}
\end{align}
$$

with $f(x) = 1.0, g(x) = 2.0$


```julia
using Gridap

# Create cartesian mesh of 10 regularly-spaced cells
L = 1.0
domain = (0.0, L, 0.0, L, 0.0, L)
n = 10
partition = (n, n, n)
model = CartesianDiscreteModel(domain, partition)

# Dirichlet BCs everywhere
labels = get_face_labeling(model)
add_tag_from_tags!(labels, "dirichlet", collect(1:8))

# Test space built from linear Lagrangian elements
order = 1
reffe = ReferenceFE(lagrangian, Float64, order)
V0    = TestFESpace(model, reffe; conformity=:H1, dirichlet_tags="dirichlet")

# Trial space
g(x) = 2.0
Ug = TrialFESpace(V0, g)

# Integration
degree = 2
Ω  = Triangulation(model)
dΩ = Measure(Ω, degree)

# Weak form
f(x)    = 1.0
h(x)    = 3.0
a(u, v) = ∫( ∇(v)⋅∇(u) )*dΩ
b(v)    = ∫( v*f )*dΩ

# Build and solve FE problem
op = AffineFEOperator(a, b, Ug, V0)
ls = LUSolver()
solver = LinearFESolver(ls)
uh = Gridap.solve(solver, op)

# Write to VTK, to be able to look at the function from Paraview
writevtk(Ω, "results", cellfields=["uh" => uh])
```

##### More details
- https://github.com/gridap/Gridap.jl
- https://gridap.github.io/Gridap.jl/dev/

## Solving PDEs using plane-wave discretisations

Consider solving the Gross-Pitaevskii equation

$$
H ψ = \left(-\frac12 Δ + V + 2 C |ψ|^2\right) ψ = μ ψ \qquad \|ψ\|_{L^2} = 1
$$

with $C = 1$ in a periodic 1D "box" of length $10$ with confining potential $V(x) = (x - 5)^2$.

This is a common mean-field model for Bose-Einstein condensates, where $\psi$ is the quantum-mechanical state of condensation and $\rho = |\psi|^2$ is the particle density.


```julia
using DFTK
using Plots

# Periodic boundary condition ... 1D problem
lattice = 10 .* [[1 0 0.];
                 [0 0 0.];
                 [0 0 0.]]

# Setup the model to be solved
V(x) = (x - 5)^2
C = 1
terms = [Kinetic(),                         # -½Δ (Kinetic energy operator in quantum physics)
         ExternalFromReal(r -> V(r[1])),
         LocalNonlinearity(ρ -> 2C * ρ^2),  # 2 * |ψ|^2
]
model = Model(lattice; n_electrons=1, terms=terms,
              spin_polarization=:spinless)  # use "spinless electrons"

# Discretise problem in a plane-wave basis
basis = PlaneWaveBasis(model; Ecut=500, kgrid=(1, 1, 1))
```


```julia
# Solve it using a constrained preconditioned LBFGS:
scfres = direct_minimization(basis, tol=1e-8);
```


```julia
# Plot the resulting particle density |ψ|^2
x = 10vec(first.(DFTK.r_vectors(basis)))
plot(x, scfres.ρ[:, 1, 1], label="ρ")
plot!(x, V.(x) / 10, label="V")
```

##### More details
- https://dftk.org
- https://docs.dftk.org

9722c6d413547f9633133404a09ed14680c776a4
ipynb
Jupyter Notebook
Projects/04_Solving_Differential_Equations.ipynb
tuanvo-git/2022-rwth-julia-workshop
0b10ede0d49b4121adb9f6932a10623f7eef7015
MIT
Projects/04_Solving_Differential_Equations.ipynb
tuanvo-git/2022-rwth-julia-workshop
0b10ede0d49b4121adb9f6932a10623f7eef7015
MIT
Projects/04_Solving_Differential_Equations.ipynb
tuanvo-git/2022-rwth-julia-workshop
0b10ede0d49b4121adb9f6932a10623f7eef7015
MIT
Qwen/Qwen-72B
1. YES
2. YES