Trying to bend my mind to Julia's way

Hi everyone,

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?

1 Like

Some ideas:

  • use enumerate to do the for loop while keeping track of the index
  • use comprehensions: M = [cos(...) for H_loco in H_range]
  • avoid global variables (for style, and for performance). Use closures to pass parameters: solve(params) = nlsolve((F, x) -> f!(F, x, params), (J, x) -> j!(J, x, params))
4 Likes

Hello @LStoleriu and welcome!

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.

https://docs.julialang.org/en/v1/manual/performance-tips/index.html (if you haven’t found it already) is good reading regarding performance.

Good luck and have fun!

4 Likes

Thank you! Good to know - this was indeed one of my concerns.

If you worry about speed (as most people seem to do here) I suggest simply trying various methods against each other with @btime from BenchmarkTools.

1 Like

Some more things:

  • 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.

Thank you - I’m looking forward to this, indeed. For now I’m just trying to be able to build different methods that work :slight_smile:

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.