# First Steps

- [DifferentialEquations.jl docs](https://diffeq.sciml.ai/dev/index.html)

## Radioactive decay

How to define your model, state variables, initial (and/or boundary) conditions, and parameters.

As a simple example, the concentration of a decaying nuclear isotope could be described as an exponential decay:

$$
\frac{d}{dt}C(t) = - \lambda C(t)
$$

**State variable(s)**
- $C(t)$: The concentration of a decaying nuclear isotope.

**Parameter(s)**
- $\lambda$: The rate constant of decay. The half-life $t_{\frac{1}{2}} = \frac{ln2}{\lambda}$

### Standard operating procedures

- Define a model function representing the right-hand-side (RHS) of the sysstem.
  - Out-of-place form: `f(u, p, t)` where `u` is the state variable(s), `p` is the parameter(s), and `t` is the independent variable (usually time). The output is the right hand side (RHS) of the differential equation system.
  - In-place form: `f!(du, u, p, t)`, where the output is saved to `du`. The rest is the same as the out of place form. The in-place form has potential performance benefits since it allocates less arrays than the out-of-place form.
- Initial conditions (`u0`) for the state variable(s).
- (Optional) parameter(s) `p`.
- Define a problem (e.g. `ODEProblem`) using the modeling function (`f`), initial conditions (`u0`), simulation time span (`tspan == (tstart, tend)`), and parameter(s) `p`.
- Solve the problem by calling `solve(prob)`.


```julia
using DifferentialEquations
using Plots

# The Exponential decay ODE model, out-of-place form
expdecay(u, p, t) = p * u

p = -1.0 # Parameter
u0 = 1.0 # Initial condition
tspan = (0.0, 2.0) # Simulation start and end time points
prob = ODEProblem(expdecay, u0, tspan, p) # Define the problem
sol = solve(prob) # Solve the problem
```

### Visualization

`DifferentialEquations.jl` has a plot recipe so that `plot(sol)` directly visualizes the time series of the state variable(s).


```julia
using Plots
plot(sol) # Visualize the solution
```

### Solution handling

Docs: https://diffeq.sciml.ai/stable/basics/solution/

`sol(t)`: solution at time `t` with interpolations.


```julia
sol(1.0)  # 
```

`sol.t`: time points of the solution. Notice *t=1* may not in one of the time points.


```julia
sol.t
```

`sol.u`: The solution at time points `sol.t`


```julia
sol.u
```

## The SIR model

A more complicated example is the [SIR model](https://www.maa.org/press/periodicals/loci/joma/the-sir-model-for-spread-of-disease-the-differential-equation-model) describing infectious disease spreading. There are more state variables and parameters.

$$
\begin{align}
\frac{d}{dt}S(t) &= - \beta S(t)I(t)  \\
\frac{d}{dt}I(t) &= \beta S(t)I(t)  - \gamma I(t)  \\
\frac{d}{dt}R(t) &= \gamma I(t)
\end{align}
$$

**State variable(s)**

- $S(t)$ : the fraction of susceptible people
- $I(t)$ : the fraction of infectious people
- $R(t)$ : the fraction of recovered (or removed) people

**Parameter(s)**

- $\beta$ : the rate of infection when susceptible and infectious people meet
- $\gamma$ : the rate of recovery of infectious people


```julia
using DifferentialEquations
using Plots

# SIR model, in-place form
function sir!(du, u, p ,t)
	s, i, r = u
	β, γ = p
	v1 = β * s * i
	v2 = γ * i
    du[1] = -v1
    du[2] = v1 - v2
    du[3] = v2
	return nothing
end

# Parameters of the SIR model
p = (β = 1.0, γ = 0.3)
u0 = [0.99, 0.01, 0.00]  # s, i, r
tspan = (0.0, 20.0)

# Define a problem
prob = ODEProblem(sir!, u0, tspan, p)

# Solve the problem
sol = solve(prob)

# Visualize the solution
plot(sol, label=["S" "I" "R"], legend=:right)
```

### Solution handling

`sol[i, j]`: `i`th component at timestep `j`


```julia
sol[2]
```


```julia
sol[1, 2]
```

`sol[i, :]`: the timeseries for the `i`th component.


```julia
sol[1, :]
```

`sol(t,idxs=1)`: the 1st element in time point(s) `t` with interpolation. `t` can be a scalar (single point) or an vector-like sequence. (multiple time points)


```julia
sol(10, idxs=2)
```


```julia
sol(0.0:0.1:20.0, idxs=2)
```

## Lorenz system

The Lorenz system is a system of ordinary differential equations having chaotic solutions for certain parameter values and initial conditions. ([Wikipedia](https://en.wikipedia.org/wiki/Lorenz_system))

$$
\begin{align}
  \frac{dx}{dt} &=& \sigma(y-x) \\
  \frac{dy}{dt} &=& x(\rho - z) -y \\
  \frac{dz}{dt} &=& xy - \beta z
\end{align}
$$

In this example, we will use [LabelledArrays.jl](https://github.com/SciML/LabelledArrays.jl) to get DSL-like syntax.


```julia
using LabelledArrays
using DifferentialEquations
using Plots

function lorenz!(du,u,p,t)
    du.x = p.σ*(u.y-u.x)
    du.y = u.x*(p.ρ-u.z) - u.y
    du.z = u.x*u.y - p.β*u.z
end

u0 = LVector(x=1.0, y=0.0, z=0.0)
p = LVector(σ=10.0, ρ=28.0, β=8/3)
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan,p)
sol = solve(prob)
```

### Visualization

`vars=(1, 2, 3)`: make a phase plot with 1st, 2nd, and the 3rd state variable. With LabelledArrays, you can use symbols instead of numbers.


```julia
plot(sol, vars=(:x, :y, :z))
```

f53773fe7edad445a73249b47d28f8ca9e107c51
ipynb
Jupyter Notebook
docs/intro/01-first-steps.ipynb
ww-juliabook/learn-diffeq
a82b9fc2d715315ec1486d94f4a8ffb98187ec9d
MIT
docs/intro/01-first-steps.ipynb
ww-juliabook/learn-diffeq
a82b9fc2d715315ec1486d94f4a8ffb98187ec9d
MIT
2022-03-19T17:02:19.000Z
2022-03-20T09:13:41.000Z
docs/01-first-steps.ipynb
ntumitolab/juliabook-diffeq
97fd4893108001d2a415291c5ccec0d2a5de54f8
MIT
Qwen/Qwen-72B
1. YES
2. YES