Created
June 14, 2014 21:59
-
-
Save acroy/28be4f2384d01f38e577 to your computer and use it in GitHub Desktop.
Demonstrates how a custom type can be used with the solvers in ODE.jl
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment