I’m a C&Maple physicist trying to learn the “true” Julia way of solving problems.
As a small project I’m trying to solve a nonlinear equation doing something like that:
using NLsolve
Nsteps = 101
Hk = 1.0
FieldMax = 3.0 * Hk
FieldMin = -FieldMax
H = FieldMax
theta_H = 45.0 * π / 180.0
theta_0 = 15.0 * π / 180.0
function f!(F, x)
F[1] = Hk * cos(x[1]) * sin(x[1]) + H * sin(x[1] - theta_H)
end
function j!(J, x)
J[1, 1] = Hk * (cos(x[1])^2 - sin(x[1])^2) + H * cos(x[1] - theta_H)
end
theta = nlsolve(f!, j!, [ 0.5 ])
println(theta.zero)
println(converged(theta))
This works great, but now I’d like to solve the same equation for Nsteps different values of H and store the results in an array called M.
I was able to do that using a C-like approach, by continuing the program like that:
M = zeros(Nsteps)
H_range = collect(range(FieldMax, length=Nsteps, stop=FieldMin))
function SWsolve(H_range, M)
i = 0
for H_loco in H_range
global H = H_loco
theta = nlsolve(f!, j!, [ 0.5 ])
i = i+1
M[i] = cos(theta_H - theta.zero[1])
end
end
SWsolve(H_range, M)
hyst = plot(H_range, M)
display(hyst)
I tried a lot to simplify this using Julia-like constructions but I failed (for example, I tried to remove the for loop from the SWsolve() function while using the dot(“.”) operator to pass the H_range elements) but I failed.
Can you suggest a more efficient way of doing this?
I’m a Julia beginner too and coming from C and python I struggle a bit with what exactly is the style of Julia. In particular, when a task can be performed by any of a for loop, broadcasting, a comprehension, or functional techniques, which is more Julian? My tentative answer is whichever feels more natural to you. Unlike python, there isn’t a performance penalty for writing your own loops and for many problems I find it easier to understand exactly what will occur using loops, and they are hard to beat for speed. The other constructions are more concise, but that doesn’t always mean faster or easier to understand.
the initial guess is always the same, so that it can be preallocated before the loop
there are functions deg2rad and rad2deg for angle conversions
as already been said, globals are bad; in some cases, though, it’s cleaner to declare some global constants, just don’t forget to add the const keyword, in which case their type is fixed, and the compiler can optimize that
Thanks!
Indeed, good point, I switched everything that’s global to const and I’ll go with deg2rad.
The initial guess is the same just for this small script posted here but as the things will get more complicated there will be some if-s and different values there.
Regarding the initial guess - unless every one must be of a different length, it’s still cheaper to pre-allocate and fill in new values on each iteration than to allocate a new array each time. In fact, allocation and GCing of a lot of short arrays may be a huge performance penalty.
Sometimes you do want to use collect, but I don’t think it’s necessary here. The range object should be perfectly fine, and in bigger cases can avoid unnecessary allocations.