March 26, 2012

simulating nonlinear oscillators using the julia language

In my post 'simulating nonlinear oscillators using octave' i showed you how to ... -- i guess you can read. Today i discovered, that a package similar to the odepkg also is included in julia. You just have to find it.

Basically it works quite similar, with some little tweaks.
There is a file called 'ode.jl' in the 'extras' folder located in the julia source folder. We just have to load this thing and have some ode solvers. The drawback is that up to now there  the only way (i found) to adapt the accuracy is to change the lines

rtol = 1.e-5
atol = 1.e-8

in the function ode23
and
tol = 1.0e-7
in the definition of oderkf{T}, which is a quite dirty method.

Another issue is that the only arguments of the ode solvers are up to now

(F::Function, tspan::AbstractVector, y_0::AbstractVector)

which means we have to hide additional parameters in our solution vector (and make its derivatives zero).

A code example, that works for me (asymmetric nonlinear potential)

function  rhs_asym(t, x)
#
# position and velocity are stored in x[1] and x[2]
#
xp = x[1]
vx = x[2]
#
# the parameters are stored in x[3] and x[4]
#
k=x[3]
alpha =x[4]
dxp = vx
dv = - (k * xp - alpha * xp^2)
dx = [dxp; dv;0;0]   # the two zeros ensure that k and alpha
# will not change
end

# problems can occur when you forget to give these values
# explicitly as Floats ...
k=1.
alpha =0.2
tend=42.

initialDisplacement = 0.0
initialVelocity     = 1.
xini = [initialDisplacement; initialVelocity; k; alpha]

#
# to get higher accuracy, change the values of rtol atol or tol
# in the ode.jl-file (or better in a copy of it).
#
# now let's integrate
(T, X) = ode45( rhs_asym, [0; tend], xini)

# we store the whole data in one matrix and save it,
# so we can load and plot with e.g. octave
#
M=[T1 Xlin]
csvwrite("oscil.csv",M)