#Understanding allocations in an RK4 calculation
1 messages · Page 1 of 1 (latest)
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
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 🌞