Why is julia updating a different struct than what I am in my function (trying to implement a 4th order runge kutta method)?

user_17022863 0
Viewed 95

In Julia I am integrating two fields in a struct: x position and x velocity. In function d(u, du) I am trying to only return the du vector without altering any values in u. u is only used to calculate the du value. Instead it is changing the values of u instead of du. I have a vector of one struct for du and one vector for u. On each time step, I update my du.x with u.xvelocity, and I update my u.xvelocity for my acceleration. For some reason it seems to break when I calculate the k2 for my runge kutta. The MWE is appended below and should run and get the same error that I am getting. Also, if my runge kutta looks incorrect, also let me know. Best.

module MyOde
mutable struct Particle 
    x::Float64
    xvelocity::Float64
end # End struct

function d(u, du) 
    for i in 1:length(u)
        testme = u[i].x 
        display(testme)
        display(u[i].xvelocity)
        du[i].x = 1.
        println("now see the issue:")
        display(u[i].x)

        testme != u[i].x ? error("\n What the heck is going on here? ") : nothing
        du[i].xvelocity = 1.
    end # End function
    return du
end # End function 


function f(u::Vector{Particle}, d, timeend, dt) 
du = Vector{Particle}(undef, length(u))
k2 = Vector{Particle}(undef, length(u))
k3 = Vector{Particle}(undef, length(u))
k4 = Vector{Particle}(undef, length(u))

for i ∈ 1:length(u)
    du[i] = Particle(0.0, 0.0)
    k2[i] = Particle(0.0, 0.0)
    k3[i] = Particle(0.0, 0.0)
    k4[i] = Particle(0.0, 0.0)
end # End list push

for i in 0.0:dt:timeend
    # Calculate the k values which will be going into the 4th order Runge-Kutta method. 
    k1 = d(u, du)
    for i ∈ 1:length(u)
        k2[i].x = u[i].x + k1[i].x *dt/2
        k2[i].xvelocity = u[i].xvelocity + k1[i].xvelocity *dt/2
    end # End k2 loop
    k2 = d(k2, du)
    for i ∈ 1:length(u)
        k3[i].x = u[i].x + k2[i].x *dt/2
        k3[i].xvelocity = u[i].xvelocity + k2[i].xvelocity *dt/2

    end # End k3 loop
    k3 = d(k3, du)
    for i ∈ 1:length(u)
        k4[i].x = u[i].x + k3[i].x *dt
        k4[i].xvelocity = u[i].xvelocity + k3[i].xvelocity *dt
    end # End k4 loop
    k4 = d(k4, du)
    for i ∈ 1:length(u)
        u[i].x += 1/6  * dt * (k1[i].x + 2k2[i].x + 2k3[i].x + k4[i].x)
        u[i].xvelocity += 1/6  * dt * (k1[i].xvelocity + 2k2[i].xvelocity + 2k3[i].xvelocity + k4[i].xvelocity)
    end # End loop 
end # End loop 
   
end # End function

u = [Particle(0.0, 0.0); Particle(6.00, 0.0)] 
timeend = .01
dt = 0.01 
@time f(u, d, timeend, dt)
end # End module
differential-equations ode julia numerical-methods struct
2 Answers

Like you said, function d(u, du) returns a mutated du. When you run k2 = d(k2, du), the variable k2 is reassigned to the same object referenced by the function variable du, and the previously referenced object is discarded, along with all the work done on it when it was still named k2 in the preceding loop.

You do the same mistake for the rest of the ks, so by the time you hit the 2nd iteration of the for i in 0.0:dt:timeend loop, the variables du, k1, k2, k3, k4 all reference the same object. The error is thrown on the 2nd iteration's k2 = d(k2, du) because in the function d, the arguments u and du reference the same object you mutate via du.

I don't know the math here so I can't give any input on that. But I'm sure you want your variables to reference separate independent objects. Look over your code and remove the assignments that cause this issue. For example, k2 = d(k2, du) should be d(k2, du).

user_41321440

It should be apparent that if you want to use du=k1,k2,k3 at the end in the update expression with their original values as slopes or derivatives, then you can not overwrite them in-between and use them as temporary points. Use one state utmp for that.

There are ways to minimize the number of arrays that need to be kept in an RK4 implementation, for instance by storing k2+k3 in k2 and using the k4 array in the role of k4, or accumulating the step update during the stages, so that only the temporary state, the current slope and the accumulated slope need to be present.

For the general problem with the flow of the contents of the data arrays see the other answer.

user_30881380
Related