7-discretizing-odes.org

discretizing odes

https://mitmath.github.io/18337/lecture7/discretizing_odes

examples

lorenz equations

function lorenz(du,u,p,t)
    du[1] = p[1]*(u[2]-u[1])
    du[2] = u[1]*(p[2]-u[3]) - u[2]
    du[3] = u[1]*u[2] - p[3]*u[3]
end
u0 = [1.0,0.0,0.0]
tspan = (0.0,20.0)
p = (10.0,28.0,8/3)
using DifferentialEquations
prob = ODEProblem(lorenz,u0,tspan,p)
sol = solve(prob)
using Plots
plot(sol)

n-body problem

using OrdinaryDiffEq
using Plots

function pleiades(du,u,p,t)
    @inbounds begin
        x = view(u,1:7)   # x
        y = view(u,8:14)  # y
        v = view(u,15:21) # x′
        w = view(u,22:28) # y′
        du[1:7] .= v
        du[8:14].= w
        for i in 15:28
            du[i] = zero(u[1])
        end
        for i=1:7,j=1:7
            if i != j
                r = ((x[i]-x[j])^2 + (y[i] - y[j])^2)^(3/2)
                du[14+i] += j*(x[j] - x[i])/r
                du[21+i] += j*(y[j] - y[i])/r
            end
        end
    end
end
tspan = (0.0,3.0)
prob = ODEProblem(pleiades,[3.0,3.0,-1.0,-3.0,2.0,-2.0,2.0,3.0,-3.0,2.0,0,0,-4.0,4.0,0,0,0,0,0,1.75,-1.5,0,0,0,-1.25,1,0,0],tspan)
sol = solve(prob,Vern8(),abstol=1e-10,reltol=1e-10)
plot(sol)

maxt = 14.0
tspan = (0.0,maxt)
prob = ODEProblem(pleiades,[3.0,3.0,-1.0,-3.0,2.0,-2.0,2.0,3.0,-3.0,2.0,0,0,-4.0,4.0,0,0,0,0,0,1.75,-1.5,0,0,0,-1.25,1,0,0],tspan)
sol = solve(prob,Vern8(),abstol=1e-10,reltol=1e-10)
plot(sol,vars=((1:7),(8:14)))

function vectorplot!(p, u, dt=.1, colors=[:orange, :blue, :green, :purple, :brown, :teal, :pink])
    x = view(u,1:7)   # x
    y = view(u,8:14)  # y
    xp = view(u,15:21) # x′
    yp = view(u,22:28) # y′

    for i in 1:7
        plot!(p, Shape([(x[i], y[i]), (x[i] + dt*xp[i], y[i] + dt*yp[i])]), linecolor=colors[i], linewidth=3)
    end

    p
end

dt = 0.1

anim = @animate for y in range(0., maxt, step=dt/2.)
    p = plot(xlim=(-10, 10), ylim=(-10, 10), markerstrokewidth=0, legend=false, foreground_color_border=:transparent,
            foreground_color_axis=:transparent, aspect_ratio=1, dpi=200, fontfamily="ETBookOT")
    #plot!(p, sol, vars=((1:7),(8:14)))
    vectorplot!(p, sol(y), dt)
end
gif(anim, "data/plieades.gif")

┌ Info: Saved animation to
│   fn = /home/radical/dev/6.338/data/plieades.gif
└ @ Plots /home/radical/.julia/packages/Plots/h3o4c/src/animation.jl:95

population ecology: lotka volterra

tspan = (0.0,200.0)
prob = ODEProblem(pleiades,[3.0,3.0,-1.0,-3.0,2.0,-2.0,2.0,3.0,-3.0,2.0,0,0,-4.0,4.0,0,0,0,0,0,1.75,-1.5,0,0,0,-1.25,1,0,0],tspan)
sol = solve(prob,Vern8(),abstol=1e-10,reltol=1e-10)
plot(sol,vars=((1:7),(8:14)))

biochemistry: robertson equations

this is actually where the classically “hard”

using Sundials, ParameterizedFunctions
function rober(du,u,p,t)
  y₁,y₂,y₃ = u
  k₁,k₂,k₃ = p
  du[1] = -k₁*y₁+k₃*y₂*y₃
  du[2] =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
  du[3] =  k₂*y₂^2
end
prob = ODEProblem(rober,[1.0,0.0,0.0],(0.0,1e5),(0.04,3e7,1e4))
sol = solve(prob,Rosenbrock23())
plot(sol)

ArgumentError: Package Sundials not found in current path:
- Run `import Pkg; Pkg.add("Sundials")` to install the Sundials package.

: :

Stacktrace:
 [1] require(::Module, ::Symbol) at ./loading.jl:876
 [2] top-level scope at In[4]:1

geometric properties

the simplest ODE is the scalar linear ODE

$$u’ = α u$$

analytic solution:

$u(t)=u(0)e^{α t}$

from this solution we have:

if $Re(α) > 0$, $lim_{t → ∞} u(t) = ∞$

if $Re(α) < 0$, $lim_{t → ∞} u(t) = 0$

if $Re(α) = 0$, solution is constant or periodic

in multivariate version $\uv’ = A\uv$

assume $A$ is diagonalizable, diagonalize:

$\uv’ = P^-1DP\uv$ $P\uv’ = DP\uv$

change coordinates $\zv = P\uv$

so we have

$\zv’ = D\zv$

but decomposed by eigenvalues. so each individual component can be treated on its own, and will behave how we talked before.

this also applies to the linearization of nonlinear ODEs.

numerically solving ODEs