Differential Equations / Apr 01 2020

Differential Equation Models

Yingbo Ma, Chris Rackauckas

If you're getting some cold feet to jump in to DiffEq land, here are some handcrafted differential equations mini problems to hold your hand along the beginning of your journey.

Classical Physics Models

Radioactive Decay of Carbon-14

using Pkgpkg"up; add DifferentialEquations"using DifferentialEquations, Plots

553.4s

Julia

First order linear ODE

$f(u, t, p) = p\frac{du}{dt}$

The Radioactive decay problem is the first order linear ODE problem of an exponential with a negative coefficient, which represents the half-life of the process in question. Should the coefficient be positive, this would represent a population growth equation.

# Half-life of Carbon-14 is 5,730 years.C₁ = 5.730# Setupu₀ = 1.0tspan = (0.0, 30.0)# Define the problemradioactivedecay(u,p,t) = -log(2) / C₁ * u#Pass to solverprob = ODEProblem(radioactivedecay, u₀, tspan)sol = solve(prob)#Plotplot(sol,linewidth=2,     title ="Carbon-14 half-life",      xaxis = "Time in thousands of years",      yaxis = "Percentage left",      label = "Numerical Solution")plot!(sol.t, t->exp(-log(2) / C₁*t), lw=3, ls=:dash,       label="Analytical Solution")

49.7s

Julia

Simple Pendulum

Second-Order Linear ODE

We will start by solving the pendulum problem. In the physics class, we often solve this problem by small angle approximation, i.e. $sin(\theta) \approx \theta$ , because otherwise, we get an elliptic integral which doesn't have an analytic solution. The linearized form is

$\ddot{\theta} + \frac{g}{L}{\theta} = 0$

But we have numerical ODE solvers! Why not solve the real pendulum?

$\ddot{\theta} + \frac{g}{L}{\sin(\theta)} = 0$

# Simple Pendulum Problemusing DifferentialEquations, Plots# Constantsconst g = 9.81const L = 1.0#Initial Conditionsu₀ = [0,π/2]tspan = (0.0,6.3)#Define the problemfunction simplependulum(du,u,p,t)    θ  = u[1]    dθ = u[2]    du[1] = dθ    du[2] = -(g/L)*sin(θ)end#Pass to solversprob = ODEProblem(simplependulum,u₀, tspan)sol = solve(prob,Tsit5())#Plotplot(sol,linewidth=2,     title ="Simple Pendulum Problem",      xaxis = "Time",      yaxis = "Height",      label = ["Theta","dTheta"])

8.1s

Julia

So now we know that behavior of the position versus time. However, it will be useful to us to look at the phase space of the pendulum, i.e., and representation of all possible states of the system in question (the pendulum) by looking at its velocity and position. Phase space analysis is ubiquitous in the analysis of dynamical systems, and thus we will provide a few facilities for it.

p = plot(sol, vars = (1,2), xlims = (-9,9),          title = "Phase Space Plot",          xaxis = "Velocity",          yaxis = "Position",          leg=false)function phase_plot(prob, u0, p, tspan=2pi)    _prob = ODEProblem(prob.f,u0,(0.0,tspan))    sol = solve(_prob, Vern9()) # Use Vern9 solver for higher accuracy    plot!(p, sol, vars = (1,2), xlims = nothing, ylims = nothing)endfor i in -4pi:pi/2:4π    for j in -4pi:pi/2:4π        phase_plot(prob, [j,i], p)    endendplot(p,xlims = (-9,9))

20.6s

Julia

Simple Harmonic Oscillator

Double Pendulum

# Double Pendulum Problemusing DifferentialEquations, Plots#Constants and setupconst m₁, m₂, L₁, L₂ = 1, 2, 1, 2initial = [0, π/3, 0, 3pi/5]tspan = (0.,50.)#Convenience function for transforming from polar to Cartesian coordinatesfunction polar2cart(sol;dt=0.02,l1=L₁,l2=L₂,vars=(2,4))    u = sol.t[1]:dt:sol.t[end]    p1 = l1*map(x->x[vars[1]], sol.(u))    p2 = l2*map(y->y[vars[2]], sol.(u))    x1 = l1*sin.(p1)    y1 = l1*-cos.(p1)    (u, (x1 + l2*sin.(p2), y1 - l2*cos.(p2)))end#Define the Problemfunction double_pendulum(xdot,x,p,t)    xdot[1]=x[2]    xdot[2]=-((g*(2*m₁+m₂)*sin(x[1])+m₂*(g*sin(x[1]-2*x[3])+2*(L₂*x[4]^2+L₁*x[2]^2*cos(x[1]-x[3]))*sin(x[1]-x[3])))/(2*L₁*(m₁+m₂-m₂*cos(x[1]-x[3])^2)))    xdot[3]=x[4]    xdot[4]=(((m₁+m₂)*(L₁*x[2]^2+g*cos(x[1]))+L₂*m₂*x[4]^2*cos(x[1]-x[3]))*sin(x[1]-x[3]))/(L₂*(m₁+m₂-m₂*cos(x[1]-x[3])^2))end#Pass to Solversdouble_pendulum_problem = ODEProblem(double_pendulum, initial, tspan)sol = solve(double_pendulum_problem, Vern7(), abs_tol=1e-10, dt=0.05);

7.1s

Julia

# Obtain coordinates in Cartesian Geometryts, ps = polar2cart(sol, l1=L₁, l2=L₂, dt=0.01)plot(ps...)

1.3s

Julia

Poincaré section

The Poincaré section is a contour plot of a higher-dimensional phase space diagram. It helps to understand the dynamic interactions and is wonderfully pretty.

The following equation came from StackOverflow question

The Poincaré section here is the collection of $(β,l_β)$ when $α=0$ and $\frac{dα}{dt}>0$ .

Hamiltonian of a double pendulum

Now we will plot the Hamiltonian of a double pendulum

using DifferentialEquations, Plots#Constants and setupinitial2 = [0.01, 0.005, 0.01, 0.01]tspan2 = (0.,200.)#Define the problemfunction double_pendulum_hamiltonian(udot,u,p,t)    α  = u[1]    lα = u[2]    β  = u[3]    lβ = u[4]    udot .=    [2(lα-(1+cos(β))lβ)/(3-cos(2β)),    -2sin(α) - sin(α+β),    2(-(1+cos(β))lα + (3+2cos(β))lβ)/(3-cos(2β)),    -sin(α+β) - 2sin(β)*(((lα-lβ)lβ)/(3-cos(2β))) + 2sin(2β)*((lα^2 - 2(1+cos(β))lα*lβ + (3+2cos(β))lβ^2)/(3-cos(2β))^2)]end# Construct a ContiunousCallbackcondition(u,t,integrator) = u[1]affect!(integrator) = nothingcb = ContinuousCallback(condition,affect!,nothing,                        save_positions = (true,false))# Construct Problempoincare = ODEProblem(double_pendulum_hamiltonian, initial2, tspan2)sol2 = solve(poincare, Vern9(), save_everystep = false, callback=cb, abstol=1e-9)function poincare_map(prob, u₀, p; callback=cb)    _prob = ODEProblem(prob.f,[0.01, 0.01, 0.01, u₀],prob.tspan)    sol = solve(_prob, Vern9(), save_everystep = false, callback=cb, abstol=1e-9)    scatter!(p, sol, vars=(3,4), markersize = 2)end

14.3s

Julia

poincare_map (generic function with 1 method)

p = scatter(sol2, vars=(3,4), leg=false, markersize = 2, ylims=(-0.01,0.03))for i in -0.01:0.00125:0.01    poincare_map(poincare, i, p)endplot(p,ylims=(-0.01,0.03))

2.8s

Julia

Hénon-Heiles System

The Hénon-Heiles potential occurs when non-linear motion of a star around a galactic center with the motion restricted to a plane.

$\begin{aligned} \frac{d^2x}{dt^2}&=-\frac{\partial V}{\partial x} \\ \frac{d^2y}{dt^2}&=-\frac{\partial V}{\partial y} \end{aligned}$

where

$V(x,y)={\frac {1}{2}}(x^{2}+y^{2})+\lambda \left(x^{2}y-{\frac {y^{3}}{3}}\right).$

We pick $\lambda=1$ in this case, so

$V(x,y) = \frac{1}{2}(x^2+y^2+2x^2y-\frac{2}{3}y^3).$

Then the total energy of the system can be expressed by

$E = T+V = V(x,y)+\frac{1}{2}(\dot{x}^2+\dot{y}^2).$

The total energy should conserve as this system evolves.

using DifferentialEquations, Plots# Setupinitial = [0.,0.1,0.5,0]tspan = (0,100.)# Remember, V is the potential of the system and T is the Total Kinetic Energy, thus E will#the total energy of the system.V(x,y) = 1//2 * (x^2 + y^2 + 2x^2*y - 2//3 * y^3)E(x,y,dx,dy) = V(x,y) + 1//2 * (dx^2 + dy^2);# Define the functionfunction Hénon_Heiles(du,u,p,t)    x  = u[1]    y  = u[2]    dx = u[3]    dy = u[4]    du[1] = dx    du[2] = dy    du[3] = -x - 2x*y    du[4] = y^2 - y -x^2end#Pass to solversprob = ODEProblem(Hénon_Heiles, initial, tspan)sol = solve(prob, Vern9(), abs_tol=1e-16, rel_tol=1e-16);

9.3s

Julia

# Plot the orbitplot(sol, vars=(1,2),      title = "The orbit of the Hénon-Heiles system",      xaxis = "x", yaxis = "y", leg=false)

1.4s

Julia

# Optional Sanity check - what do you think this returns and why?println(sol.retcode)

0.3s

Julia

#Plot -plot(sol, vars=(1,3),      title = "Phase space for the Hénon-Heiles system",      xaxis = "Position", yaxis = "Velocity")plot!(sol, vars=(2,4), leg = false)

1.1s

Julia

# We map the Total energies during the time intervals of the solution (sol.u here) to a new vector# pass it to the plotter a bit more convenientlyenergy = map(x->E(x...), sol.u)ΔE = energy[1]-energy[end]

0.4s

Julia

-3.0986e-5

# Plotplot(sol.t, energy,      title = "Change in Energy over Time",      xaxis = "Time in iterations",      yaxis = "Change in Energy")

0.3s

Julia

Symplectic Integration

To prevent energy drift, we can instead use a symplectic integrator. We can directly define and solve the SecondOrderODEProblem:

function HH_acceleration!(dv,v,u,p,t)    x,y  = u    dx,dy = dv    dv[1] = -x - 2x*y    dv[2] = y^2 - y -x^2endinitial_positions = [0.0,0.1]initial_velocities = [0.5,0.0]prob = SecondOrderODEProblem(           HH_acceleration!,           initial_velocities,           initial_positions,           tspan)sol2 = solve(prob, KahanLi8(), dt=1/10);

7.9s

Julia

Notice that we get the same results:

# Plot the orbitplot(sol2, vars=(3,4),      title = "The orbit of the Hénon-Heiles system",      xaxis = "x", yaxis = "y", leg=false)

2.2s

Julia

plot(sol2, vars=(3,1),      title = "Phase space for the Hénon-Heiles system",      xaxis = "Position", yaxis = "Velocity")plot!(sol2, vars=(4,2), leg = false)

0.9s

Julia

but now the energy change is essentially zero:

energy = map(x->E(x[3], x[4], x[1], x[2]), sol2.u)ΔE = energy[1]-energy[end]

0.3s

Julia

9.02056e-15

# GKS: Rectangle definition is invalid in routine SET_WINDOW# Since energy is essentially zeroplot(sol2.t, energy, title = "Change in Energy over Time",     xaxis = "Time in iterations", yaxis = "Change in Energy")

0.5s

Julia

It's so close to zero it breaks GR! And let's try to use a Runge-Kutta-Nyström solver to solve this. Note that Runge-Kutta-Nyström isn't symplectic.

sol3 = solve(prob, DPRKN6());energy = map(x->E(x[3], x[4], x[1], x[2]), sol3.u)ΔE = energy[1]-energy[end]

5.3s

Julia

-8.83687e-6

plot(sol3.t, energy, title = "Change in Energy over Time",      xaxis = "Time in iterations", yaxis = "Change in Energy")

0.3s

Julia

Note that we are using the DPRKN6 solver at reltol=1e-3 (the default), yet it has a smaller energy variation than Vern9 at abs_tol=1e-16, rel_tol=1e-16. Therefore, using specialized solvers to solve its particular problem is very efficient.

The Outer Solar System

DiffEqPhysics
RecursiveArrayTools

pkg"up; add DiffEqPhysics RecursiveArrayTools"using Plots, DifferentialEquations, DiffEqPhysics, RecursiveArrayTools

32.6s

Julia

G = 2.95912208286e-4M = [1.00000597682, 0.000954786104043, 0.000285583733151, 0.0000437273164546, 0.0000517759138449, 1/1.3e8]planets = ["Sun", "Jupiter", "Saturn", "Uranus", "Neptune", "Pluto"]pos_x = [0.0,-3.5023653,9.0755314,8.3101420,11.4707666,-15.5387357]pos_y = [0.0,-3.8169847,-3.0458353,-16.2901086,-25.7294829,-25.2225594]pos_z = [0.0,-1.5507963,-1.6483708,-7.2521278,-10.8169456,-3.1902382]pos = ArrayPartition(pos_x,pos_y,pos_z)vel_x = [0.0,0.00565429,0.00168318,0.00354178,0.00288930,0.00276725]vel_y = [0.0,-0.00412490,0.00483525,0.00137102,0.00114527,-0.00170702]vel_z = [0.0,-0.00190589,0.00192462,0.00055029,0.00039677,-0.00136504]vel = ArrayPartition(vel_x,vel_y,vel_z)tspan = (0.,200_000)

0.6s

Julia

(0.0, 200000)

The N-body problem's Hamiltonian is

Here, we want to solve for the motion of the five outer planets relative to the sun, namely, Jupiter, Saturn, Uranus, Neptune and Pluto.

const ∑ = sumconst N = 6potential(p, t, x, y, z, M) = -G*∑(i->∑(j->(M[i]*M[j])/sqrt((x[i]-x[j])^2 + (y[i]-y[j])^2 + (z[i]-z[j])^2), 1:i-1), 2:N)

0.4s

Julia

potential (generic function with 1 method)

Hamiltonian System

NBodyProblem constructs a second order ODE problem under the hood. We know that a Hamiltonian system has the form of

For an N-body system, we can symplify this as:

Thus $\dot{q}$ is defined by the masses. We only need to define $\dot{p}$ , and this is done internally by taking the gradient of V. Therefore, we only need to pass the potential function and the rest is taken care of.

nprob = NBodyProblem(potential, M, pos, vel, tspan)sol_nbody = solve(nprob, Yoshida6(), dt=100)

10.8s

Julia

# Weird resultsorbitplot(sol_nbody, body_names=planets)

3.7s

Julia

Kepler Problem and Sympletic solvers

The Hamiltonian $\mathcal {H}$ and the angular momentum $L$ for the Kepler problem are

Also, we know that

using Pkgpkg"up; add ForwardDiff"using DifferentialEquations, LinearAlgebra, ForwardDiff, Plots

5.1s

Julia

H(q,p) = norm(p)^2/2 - inv(norm(q))AM(q,p) = q[1]*p[2] - p[1]*q[2]pdot(dp,p,q,params,t) = ForwardDiff.gradient!(dp, q->-H(q, p), q)qdot(dq,p,q,params,t) = ForwardDiff.gradient!(dq, p-> H(q, p), p)initial_position = [.4, 0]initial_velocity = [0., 2.]initial_cond = (initial_position, initial_velocity)initial_first_integrals = (H(initial_cond...), AM(initial_cond...))tspan = (0,20.)prob = DynamicalODEProblem(pdot, qdot, initial_velocity, initial_position, tspan)sol = solve(prob, KahanLi6(), dt=1//10);

5.3s

Julia

Let's plot the orbit and check the energy and angular momentum variation. We know that energy and angular momentum should be constant, and they are also called first integrals.

plot_orbit(sol) = plot(sol,vars=(3,4), lab="Orbit", title="Kepler Problem Solution")function plot_first_integrals(sol, H, AM)    plot(initial_first_integrals[1].-map(u->H(u[2,:], u[1,:]), sol.u), lab="Energy variation", title="First Integrals")    plot!(initial_first_integrals[2].-map(u->AM(u[2,:], u[1,:]), sol.u), lab="Angular momentum variation")endanalysis_plot(sol, H, AM) = plot(plot_orbit(sol), plot_first_integrals(sol, H, AM))

0.8s

Julia

analysis_plot (generic function with 1 method)

analysis_plot(sol, H, AM)

4.4s

Julia

Let's try to use a Runge-Kutta-Nyström solver to solve this problem and check the first integrals' variation.

sol2 = solve(prob, DPRKN6())  # dt is not necessary, because unlike symplectic                              # integrators DPRKN6 is adaptivelength(sol2.u)

2.9s

Julia

analysis_plot(sol2, H, AM)

1.7s

Julia

We can see that ERKN4 does a bad job for this problem, because this problem is not sinusoid-like.

One advantage of using DynamicalODEProblem is that it can implicitly convert the second-order ODE problem to a normal system of first-order ODEs, which is solvable for other ODE solvers. Let's use the Tsit5 solver for the next example.

sol4 = solve(prob, Tsit5())length(sol4.u)

8.1s

Julia

analysis_plot(sol4, H, AM)

1.7s

Julia

Note: There is drifting for all the solutions, and high order methods are drifting less because they are more accurate.

Conclusion

Symplectic integrator does not conserve energy completely at all times, but the energy can come back. In order to make sure that the energy fluctuation comes back eventually, the symplectic integrator has to have a fixed time step. Despite the energy variation, the symplectic integrator conserves the angular momentum perfectly.

Both Runge-Kutta-Nyström and Runge-Kutta integrator do not conserve energy nor the angular momentum, and the first integrals do not tend to come back. An advantage Runge-Kutta-Nyström integrator over symplectic integrator is that RKN integrator can have adaptivity. An advantage Runge-Kutta-Nyström integrator over Runge-Kutta integrator is that RKN integrator has less function evaluation per step. The ERKN4 solver works best for sinusoid-like solutions.

Manifold Projection

In this example, we know that energy and angular momentum should be conserved. We can achieve this through manifold projection. As the name implies, it is a procedure to project the ODE solution to a manifold. Let's start with a base case, where manifold projection isn't being used.

using Pkgpkg"up; add DiffEqCallbacks"using DiffEqCallbacks, DifferentialEquations

5.1s

Julia

plot_orbit2(sol) = plot(sol,vars=(1,2), lab="Orbit", title="Kepler Problem Solution")function plot_first_integrals2(sol, H, AM)    plot(initial_first_integrals[1].-map(u->H(u[1:2],u[3:4]), sol.u), lab="Energy variation", title="First Integrals")    plot!(initial_first_integrals[2].-map(u->AM(u[1:2],u[3:4]), sol.u), lab="Angular momentum variation")endanalysis_plot2(sol, H, AM) = plot(plot_orbit2(sol), plot_first_integrals2(sol, H, AM))function hamiltonian(du,u,params,t)    q, p = u[1:2], u[3:4]    qdot(@view(du[1:2]), p, q, params, t)    pdot(@view(du[3:4]), p, q, params, t)endprob2 = ODEProblem(hamiltonian, [initial_position; initial_velocity], tspan)sol_ = solve(prob2, RK4(), dt=1//5, adaptive=false)analysis_plot2(sol_, H, AM)

6.4s

Julia

There is a significant fluctuation in the first integrals, when there is no manifold projection.

function first_integrals_manifold(residual,u)    residual[1:2] .= initial_first_integrals[1] - H(u[1:2], u[3:4])    residual[3:4] .= initial_first_integrals[2] - AM(u[1:2], u[3:4])endcb = ManifoldProjection(first_integrals_manifold)sol5 = solve(prob2, RK4(), dt=1//5, adaptive=false, callback=cb)analysis_plot2(sol5, H, AM)

16.0s

Julia

We can see that thanks to the manifold projection, the first integrals' variation is very small, although we are using RK4 which is not symplectic. But wait, what if we only project to the energy conservation manifold?

function energy_manifold(residual,u)    residual[1:2] .= initial_first_integrals[1] - H(u[1:2], u[3:4])    residual[3:4] .= 0endenergy_cb = ManifoldProjection(energy_manifold)sol6 = solve(prob2, RK4(), dt=1//5, adaptive=false, callback=energy_cb)analysis_plot2(sol6, H, AM)

4.1s

Julia

There is almost no energy variation but angular momentum varies quite bit. How about only project to the angular momentum conservation manifold?

function angular_manifold(residual,u)    residual[1:2] .= initial_first_integrals[2] - AM(u[1:2], u[3:4])    residual[3:4] .= 0endangular_cb = ManifoldProjection(angular_manifold)sol7 = solve(prob2, RK4(), dt=1//5, adaptive=false, callback=angular_cb)analysis_plot2(sol7, H, AM)

3.3s

Julia

Again, we see what we expect.

Bayesian Inference on a Pendulum using Turing.jl

Vaibhav Dixit

Let's define our simple pendulum problem. Here our pendulum has a drag term ω and a length L. We get first-order equations by defining the first term as the velocity and the second term as the position, getting

using Pkgpkg"up; add DiffEqBayes RecursiveArrayTools StatsPlots Distributions Turing@0.10.1"using DiffEqBayes, DifferentialEquations, RecursiveArrayTools, Distributions, Plots, StatsPlots

97.6s

Julia

function pendulum(du,u,p,t)    ω,L = p    x,y = u    du[1] = y    du[2] = - ω*y -(9.8/L)*sin(x)endu0 = [0.0,0.1]tspan = (0.0,10.0)prob1 = ODEProblem(pendulum,u0,tspan,[1.0,2.5])

315.2s

Julia

ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true timespan: (0.0, 10.0) u0: [0.0, 0.1]

Solve the model and plot

To understand the model and generate data, let's solve and visualize the solution with the known parameters

sol = solve(prob1,Tsit5())plot(sol)

22.1s

Julia

It's the pendulum, so you know what it looks like. It's periodic, but since we have not made a small angle assumption it's not exactly sin or cos. Because the true dampening parameter ω is 0, the solution does not decay over time, nor does it increase. The length L determines the period.

Create some dummy data to use for estimation

We now generate some dummy data to use for estimation

t = collect(range(1,stop=10,length=10))randomized = VectorOfArray([(sol(t[i]) + 0.01randn(2)) for i in 1:length(t)])data = convert(Array,randomized)

1.6s

Julia

2×10 Array{Float64,2}: 0.0242092 -0.000847219 -0.00756813 … 0.00131001 0.000578024 -0.0222702 -0.0271954 0.0409507 0.00534257 0.00919522

Let's see what our data looks like on top of the real solution

scatter!(data')

3.3s

Julia

This data captures the non-dampening effect and the true period, making it perfect to attempting a Bayesian inference.

Perform Bayesian Estimation

Now let's fit the pendulum to the data. Since we know our model is correct, this should give us back the parameters that we used to generate the data! Define priors on our parameters. In this case, let's assume we don't have much information but have a prior belief that ω is between 0.1 and 3.0, while the length of the pendulum L is probably around 3.0:

priors = [Uniform(0.1,3.0), Normal(3.0,1.0)]

1.5s

Julia

2-element Array{Distribution{Univariate,Continuous},1}: Uniform{Float64}(a=0.1, b=3.0) Normal{Float64}(μ=3.0, σ=1.0)

Finally, let's run the estimation routine from DiffEqBayes.jl using the Turing.jl backend.

using DiffEqBayesbayesian_result = turing_inference(prob1,Tsit5(),t,data,priors;                       num_samples=10_000, syms = [:ω,:L])

68.4s

Julia

Notice that while our guesses had the wrong means, the learned parameters converged to the correct means, meaning that it learned good posterior distributions for the parameters. To look at these posterior distributions on the parameters, we can examine the chains:

plot(bayesian_result)

11.6s

Julia

As a diagnostic, we will also check the parameter chains. The chain is the MCMC sampling process. The chain should explore parameter space and converge reasonably well, and we should be taking a lot of samples after it converges (it is these samples that form the posterior distribution!)

plot(bayesian_result, colordim = :parameter)

1.8s

Julia

Notice that after a while these chains converge to a "fuzzy line", meaning it found the area with the most likely and then starts to sample around there, which builds a posterior distribution around the true mean.

Conditional Dosing Pharmacometric Example

using DifferentialEquations, Plots, DiffEqCallbacks

Shift+Enter to run

Julia

In this example, we will show how to model conditional dosing using the DiscreteCallbacks. The problem is as follows. The patient has a drug A(t) in their system. The concentration of the drug is given as C(t)=A(t)/V for some volume constant V. At t=4, the patient goes to the clinic and is checked. If the concentration of the drug in their body is below 4, then they will receive a new dose.

For our model, we will use the simple decay equation. We will write this in the in-place form to make it easy to extend to more complicated examples

function f!(du,u,p,t)    du[1] = p * u[1]endu0 = [10.0]const V2 = 1prob = ODEProblem(f!,u0,(0.0,10.0), -1.0)

0.8s

Julia

ODEProblem with uType Array{Float64,1} and tType Float64. In-place: true timespan: (0.0, 10.0) u0: [10.0]

Let's see what the solution looks like without any events.

sol = solve(prob,Tsit5())plot(sol)

3.8s

Julia

We see that at time t=4, the patient should receive a dose. Let's code up that event. We need to check at t=4 if the concentration u[1]/4 is <4, and if so, add 10 to u[1]. We do this with the following:

condition(u,t,integrator) = t==4 && u[1]/V2<4affect!(integrator) = integrator.u[1] += 10cb = DiscreteCallback(condition, affect!)

0.8s

Julia

DiscreteCallback{typeof(condition),typeof(affect!),typeof(INITIALIZE_DEFAULT)}(condition, affect!, INITIALIZE_DEFAULT, Bool[1, 1])

Now we will give this callback to the solver, and tell it to stop at t=4 so that way the condition can be checked

sol = solve(prob,Tsit5(),tstops=[4.0],callback=cb)plot(sol)

2.5s

Julia

Let's show that it actually added 10 instead of setting the value to 10. We could have set the value using affect!(integrator) = integrator.u[1] = 10

println(sol(4.00000))println(sol(4.000000000001))

0.9s

Julia

Now let's model a patient whose decay rate for the drug is lower:

u0 = [10.0]prob = ODEProblem(f!,u0,(0.0,10.0), -1/6)sol = solve(prob,Tsit5())plot(sol)

0.9s

Julia

Under the same criteria, with the same event, this patient will not receive a second dose:

sol = solve(prob,Tsit5(),tstops=[4.0],callback=cb)plot(sol)

0.6s

Julia

Biochemical reaction model

Please see DiffEqBiological.jl examples.

Differential Equation Models

Classical Physics Models

Radioactive Decay of Carbon-14

Simple Pendulum

Simple Harmonic Oscillator

Poincaré section

Hénon-Heiles System

The Outer Solar System

Kepler Problem and Sympletic solvers

Manifold Projection

Bayesian Inference on a Pendulum using Turing.jl

Solve the model and plot

Create some dummy data to use for estimation

Perform Bayesian Estimation

Conditional Dosing Pharmacometric Example

Biochemical reaction model

Runtimes (1)