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)

# add the path to your juliasource/extras folder

load("/Users/user/juliadir/extras/ode.jl")

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)

## No comments:

## Post a Comment