My code solves the Monte Carlo problem of a large ensemble of particles interacting with field. I use the Newton’s Equations to evolve the system. Numerically, ordinary equation is needed. Some pieces of my code are as follows:
function Eq1(t, y) (x, y, px, py) = y r = sqrt( x^2 + y^2 + s1 ) x_prime = px y_prime = py px_prime = -2x/r^3 py_prime = -2y/r^3 [x_prime; y_prime; px_prime; py_prime] end function Init2() x = Array(Float64, ne) y = similar(x) px = similar(x) py = similar(x) rand1 = rand(ne) x .= R_init * sin(2pi.*rand1) y .= R_init * cos(2pi.*rand1) rand1 = rand(ne) px .= p_init * sin(2pi.*rand1) py .= p_init * cos(2pi.*rand1) println(R_init) println(p_init) println(s1) for i in 1:ne tout, yout = ode45(Eq1, [x[i], y[i], px[i], py[i]], [0.0; 100.0])#; points =:specified) x[i] = yout y[i] = yout px[i] = yout py[i] = yout if mod(i, 1000) == 0 println("i=$(i)") end end
Whe I set the parameters as
R0_init=0.5 p0_init=0.7869 s1 = 0.5 ne = 5000
the @time Inite() command give a result as:
14.935459 seconds (267.48 M allocations: 8.997 GB, 4.63% gc time)
I have solved similar problem with Fortran and Matlab, both of them are much more efficient. I guess I have wrongly used the package, but I can only find very little details in the package github page. Can anyone tell me what the problem is and how can i fix it?