Mudanças na v3 do OrdinaryDiffEq.jl

Olá!

Eu sei, um longo, longo tempo passou desde a minha última postagem. Eu gostaria de me desculpar! Estive muito ocupado devido a uma quantidade enorme de trabalho a ser feito no INPE.

De qualquer forma, esse será uma postagem pequena. Apenas gostaria de atualizar o código-fonte da última lição de tal forma que rode corretamente utilizando Julia v0.6 e OrdinaryDiffEq.jl v3.

Como eu disse, o pacote OrdinaryDiffEq.jl, que é parte do projeto JuliaDiffEq, é maravilhoso. Na minha humilde opinião, é um dos melhores solvers de equações diferenciais de código aberto disponíveis. O desenvolvimento dessa suíte está ocorrendo a uma taxa muito elevada. Isso é bom, pois o pacote tem ficado melhor a cada dia. O lado ruim é que a API também se altera muito rápido e às vezes seu código precisa ser modificado para continuar a funcionar nas novas versões.

Experiência necessária

Isso é apenas uma atualização do código-fonte apresentado na minha última postagem. Portanto, recomendo fortemente a lê-la primeiro:

Mudanças na v3 do OrdinaryDiffEq.jl

Desde que eu postei meu último tutorial, a API do OrdinaryDiffEq.jl mudou (na verdade, a API é definida em DiffEqBase.jl, mas ele será atualizado quando você atualizar o OrdinaryDiffEq.jl). Existem três modificações que devem ser feitas para que o código-fonte anterior possa funcionar corretamente com Julia v0.6 e OrdinaryDiffEq.jl v3.0.3.

Interface de Array de Dados Customizada

Anteriormente, uma interface de array de dados customizada podia ser criada através de uma nova estrutura que herdava o tipo abstrato DEDataArray{T}. Isso se alterou e agora o novo tipo abstrato correto é DEDataVector{T}. Então, nossa definição do SimWorkspace deve ser:

type SimWorkspace{T} <: DEDataVector{T}
    x::Array{T,1}
    a::T
end

Assinatura da equação dinâmica

Na versão utilizada no último tutorial, a assinatura da função dinâmica era:

function f(t,u,du)

Devido a algumas mudanças internas, a assinatura agora é:

function f(du,u,t,p)

p é um vetor de parâmetros que não é utilizado nesse exemplo, mas simplifica muitas coisas. Se você quiser saber mais informações, por favor veja http://docs.juliadiffeq.org/stable/tutorials/ode_example.html#Defining-Parameterized-Functions-1.

Assinatura da função de condição do callback

A assinatura da função de condição do callback se alterou de:

function condition_control_loop(t,u,integrator)

para:

function condition_control_loop(u,t,integrator)

Código-fonte atualizado

Finalmente, o código-fonte atualizado do último tutorial pode ser visto a seguir. Ele irá fornecer exatamente o mesmo resultado.

using DiffEqBase
using OrdinaryDiffEq

###############################################################################
#                                  Variables
###############################################################################

# Global variable to store the acceleration.
# This solution was selected just for the sake of simplification. Don't use
# global variables!
a = 0.0

# Parameters.
r = 10.0   # Reference position.
k = 0.3    # Proportional gain.
d = 0.8    # Derivative gain.

# Configuration of the simulation.
tf     = 30.0             # Final simulation time.
tstops = collect(0:1:tf)  # Instants that the control loop will be computed.
u0     = [0.0; 0.0]       # Initial state.

###############################################################################
#                                  Functions
###############################################################################

# Dynamic equation.
function dyn_eq(du,u,p,t)
    du .= [0 1; 0 0]*u + [0;1]*a
end

# CALLBACK: Control loop.
# =======================

# Condition function.
function condition_control_loop(u,t,integrator)
    (t in tstops)
end

# Affect function.
function control_loop!(integrator)
    global a

    p = integrator.u[1]
    v = integrator.u[2]

    a = k*(r-p) + d*(0.0-v)
end

cb   = DiscreteCallback(condition_control_loop, control_loop!)

###############################################################################
#                                    Solver
###############################################################################

prob = ODEProblem(dyn_eq, u0, (0.0, tf))
sol  = solve(prob, Tsit5(), callback = cb, tstops=tstops)
nothing

Conclusões

Esse tutorial apenas atualiza o código-fonte do último de tal forma que se possa executá-lo utilizando Julia v0.6 e OrdinaryDiffEq.jl v3.0.3.

Estou planejando para o meu próxima tutorial um exemplo simples de como codificar filtros de Kalman em linguagem Julia!