#Understanding allocations in an RK4 calculation

1 messages · Page 1 of 1 (latest)

sand vortex
#

I think
ypts[i,:] is to blame
I think it allocates a new array each time (so, at least 5 times/point)
Try using view function
🌞

EDIT: also, I think my_func_2o allocates too
Try mutating y instead of allocating a new array

#

something like

dummy_variable = y[1]
y[1] = y[2]
y[2] = -dummy_variable
return y

#

So, you don't make a new array every time you call your function

Using "view" in the solver should reduce allocation even more

#

I don't have access to my PC atm, so I can't test it myself

But changing ypts[i,:] to view(ypts,i,:) might do the trick

EDIT: Ok, I think it might not solve all problems...
I'll be able to look into it when I'm back home 🫡
Or perhaps, someone else can help before that

sand vortex
#

Sorry, I shouldn't have answered without properly waking up

Here's a version without allocation

function rk4_2o(f, init::Vector{Float64}, Npts::Int64, t::Float64)
    n = length(init)
    tpts = LinRange(0,t,Npts)
    h = tpts[2]
    ypts = Matrix{Float64}(undef,Npts,n)
    ypts[1, :] = init

    placeholder = Vector{Float64}(undef,n)
    k1 = Vector{Float64}(undef,n)
    k2 = Vector{Float64}(undef,n)
    k3 = Vector{Float64}(undef,n)
    k4 = Vector{Float64}(undef,n)

    for i in 1:(Npts-1)
        f(k1,tpts[i],view(ypts,i,:))
        for j in 1:n
            placeholder[j] = ypts[i,j] + h*k1[j]/2
        end
        f(k2,tpts[i]+0.5h,placeholder)
        for j in 1:n
            placeholder[j] = ypts[i,j] + h*k2[j]/2
        end
        f(k3,tpts[i] + 0.5*h,placeholder)
        for j in 1:n
            placeholder[j] = ypts[i,j] + h*k3[j]
        end
        f(k4,tpts[i] + h, placeholder)
        for j in 1:n
            ypts[i+1,j] = ypts[i,j] + h*(k1[j]+2k2[j]+2k3[j]+k4[j])/6
        end
    end
    return tpts, ypts
end
function my_func_2o!(dy,t, y)
    dy[1],dy[2] = y[2],-y[1]
    return dy
end

y = [1.0,0.0]
t,sol = rk4_2o(my_func_2o!,y,1000,10.0)
#

// I mean... 6 allocations, just like in your 1st order version

So, instead of creating new arrays every time, there's a "placeholder" array that saves function evaluation results
that's why we added a new argument to "my_func_2o" ... so we don't modify ypts but just calculated things based off ypts

You can get rid of these for loop inside 🌞