Demonstrates how a custom type can be used with the solvers in ODE.jl

ode_custom_type.jl

using ODE

import Base.norm  # we will extend norm to support our new type

const delta0 = 0.
const V0 = 1.
const g0 = 0.

################################################################################
# define custom type ...
type CompSol
  rho::Matrix{Complex128}
  x::Float64
  p::Float64

  CompSol(r,x,p) = new(copy(r),x,p)
end

# ... which has to support the following operations
# to work with odeX
norm(y::CompSol, p::Float64) = maximum([Base.norm(y.rho, p) abs(y.x) abs(y.p)])

+(y1::CompSol, y2::CompSol) = CompSol(y1.rho+y2.rho, y1.x+y2.x, y1.p+y2.p)
-(y1::CompSol, y2::CompSol) = CompSol(y1.rho-y2.rho, y1.x-y2.x, y1.p-y2.p)
*(y1::CompSol, s::Real) = CompSol(y1.rho*s, y1.x*s, y1.p*s)
*(s::Real, y1::CompSol) = y1*s
/(y1::CompSol, s::Real) = CompSol(y1.rho/s, y1.x/s, y1.p/s)
################################################################################

# define RHSs of differential equations
# delta, V and g are parameters
function rhs(t, y, delta, V, g)
  H = [[-delta/2 V]; [V delta/2]]

  rho_dot = -im*H*y.rho + im*y.rho*H
  x_dot = y.p
  p_dot = -y.x

  return CompSol( rho_dot, x_dot, p_dot)
end

# inital conditons
rho0 = zeros(2,2);
rho0[1,1]=1.;
y0 = CompSol(complex(rho0), 2., 1.)

# solve ODEs
t,y = ODE.ode45((t,y)->rhs(t, y, delta0, V0, g0), y0, [0., 2pi]);


using Winston
# dynamics of rho
plot(t,  Float64[real(R.rho[1,1]) for R in y], t, cos(t).^2, "bo")
oplot(t, Float64[real(R.rho[2,2]) for R in y], "--", t, sin(t).^2, "rs")

# dynamics of x and p
plot(t,  Float64[R.x for R in y], t, 2*cos(t)+sin(t), "bo")
oplot(t, Float64[R.p for R in y], "--", t, cos(t)-2*sin(t), "rs")