[julia-users] Re: Ode solver thought, numerical recipes

[François Fayard]

Hi,

Here is a first shot at the program (My first Julia program). It seems that 
I don't get good performance (46% of time spent in gc). What is the reason 
for that?

type EulerSolver
  f::Function
  t0::Float64
  y0::Vector{Float64}
  delta_t::Float64
  delta_t_save::Float64
  t_right::Vector{Float64}
  y_right::Vector{Vector{Float64}}
  t_left::Vector{Float64}
  y_left::Vector{Vector{Float64}}
end

function ode_solver(f, t0::Float64, y0::Vector{Float64}, delta_t::Float64, 
delta_t_save::Float64)
  solver = EulerSolver(f, t0, y0, delta_t, delta_t_save, [], [], [], [])
  return solver
end

function evaluate(solver::EulerSolver, t_final::Float64)
  t::Float64
  y::Vector{Float64}
  if length(solver.t_right) > 0
    if t_final >= solver.t_right[length(solver.t_right)]
      t = solver.t_right[length(solver.t_right)]
      y = solver.y_right[length(solver.y_right)]
    elseif t_final < solver.t_right[1]
      t = solver.t0
      y = solver.y0
    else
      k = 1
      while t >= solver.t_right[k]
        k += 1
      end
      t = solver.t_right[k]
      y = solver.y_right[k]
    end
  else
    t = solver.t0
    y = solver.y0
  end
  
  last_t_saved::Float64 = t
  delta_t::Float64 = solver.delta_t
  dy_dt::Vector{Float64} = zeros(y)

  is_finished::Bool = false
  while !is_finished
    if t + delta_t > t_final
      delta_t = t_final - t
      is_finished = true
    end

    f!(dy_dt, t, y)
    for k = 1:length(y)
      y[k] = y[k] + delta_t * dy_dt[k]
    end
    t = t + delta_t

    if abs(t - last_t_saved) >= solver.delta_t_save
      println("Save")
      last_t_saved = t
      push!(solver.t_right, t)
      push!(solver.y_right, deepcopy(y))
    end
  end

  return y
end

function f!(dy_dt::Vector{Float64}, t::Float64, y::Vector{Float64})
  dy_dt[1] = y[1]
end

t0 = 0.0
y0 = [ 1.0 ]
delta_t = 1.0e-8
delta_t_save = 1.0

solver = ode_solver(f!, t0, y0, delta_t, delta_t_save)
@time println(evaluate(solver, 3.12))