Getting Stack Overflow error when running Derivative(x)^n

I have published my code in the following link.

I am pasting my code here as well

##
using Symbolics
using Plots
using Statistics
##
a = 0.529*10^-10 # Radius of hydrogen atom in meters
function Legrande(l::Int, cosθ::Num)
    # Compute the Legendre polynomial. Here cosθ is a symbolic variable
    D_cst_l = Differential(cosθ)^l # This is lᵗʰ derivative operator wrt cst
    legrande = (1 / (2^l * factorial(l))) * (D_cst_l(cosθ^2 - 1)^l)
    legrande = expand_derivatives(legrande)
    return legrande
end
#%%
function AssoLegrande(l::Int, m::Int, cst::Float64)
    # Define the symbolic variable
    @variables cosθ
    # Compute the associated Legendre polynomial
    D_cst_m = Differential(cosθ)^m
    asslegrande = (-1)^m * (1 - cosθ^2)^(m/2) * D_cst_m(Legrande(l, cosθ))
    asslegrande = expand_derivatives(asslegrande)
    # Substitute the numeric value of cst into the symbolic expression
    asslegrande = substitute(asslegrande, cosθ => cst)
    return asslegrande
end
#%%
function Laguerre(q::Int,x::Num)
    # Laguerre function
    # Here x is a symbolic variable
    if q==0
        leg=1;
    else
        Dx_q = Differential(x)^q
        leg = exp(x)*Dx_q(exp(-x)*x^q)/factorial(q)
        leg = expand_derivatives(leg)
    end
    return leg
    # this returns Laguerre polynomial
end
#%%
function AssoLaguerre(p ::Int, q ::Int, k::Float64)
    # Associated Laguerre function
    @variables x
    if p == 0
        Dx = Differential(x)
        asso_lag = (-1)^p * (p+q,x)
    else
        Dx_p = Differential(x)^p
        asso_lag = (-1)^p * Dx_p(Laguerre(p+q, x))
    end
    # this returns Associated Laguerre polynomial
    asso_lag = expand_derivatives(asso_lag)
    asso_lag = substitute(asso_lag, x => k)
    return asso_lag
end
#%%
function Radial(n::Int, l::Int, r::Float64)
    # Define the symbolic variables
    @variables ρ
    # Compute the radial part of the wavefunction
    radial = sqrt((2/(n*a))^3 * factorial(n-l-1) / (factorial(n+l) * 2 * n)) * exp(-r/(n*a)) * (2*r/(n*a))^l
    radial = radial * AssoLaguerre(2*l+1, n-l-1, r)
    radial = expand_derivatives(radial)
    # Substitute temp = 2*r/(n*a)
    ρ_val = 2*r/(n*a)
    
    # Substitute the numeric value of r into the symbolic expression
    radial = substitute(radial,Dict([ρ => ρ_val]))
    
    return radial
end
#%%
function Angular(l::Int, m::Int, cst::Float64, phi::Float64)
    i = Complex{Float64}(0, 1)  # Define the imaginary unit
    angular = sqrt((2*l+1) * factorial(l-m) / (4*pi * factorial(l+m))) * AssoLegrande(l, m, cst) * exp(i * m * phi)
    return angular
end
##
@variables ρ, cosθ, φ # ρ is radial distance, φ is azimuthal angle, θ is the 
# n, l, m are the principle, azimuthal, and magnetic quantum numbers
##% Initialisations
rₘᵢₙ = 0.5* a
rₘₐₓ = 50* a
r_steps = 10
θ_steps = 10
φ_steps = 10
Total_steps = r_steps * θ_steps * φ_steps
n = 1
l = 0
m = 0
x = zeros(Total_steps)
y = x
z = x
#ψ(n,l,m,θ,r,ϕ) = Radial(n, l, r)*Angular(l, m, cos(θ), ϕ)
# ψ = Radial(...) * Angular(...)
# The Radial(...) calls the Associated Laguerre function which in turn calls the Laguerre function
# Simillarly the Radial(...) calls the Associated Legendre function which in turn calls the Legendre function
ψ² = zeros(Total_steps) # Initializing place holder for |ψ|²
global  counter = 1
##
for r = range(rₘᵢₙ, rₘₐₓ, 100)
    for θ = range(0, π, 30)
        for ϕ=range(0, 2*π, 30)
            print(ϕ)
            # substitute and save the values of psi for all
            # different psi, theta and r
            ψ = Radial(n, l, r)*Angular(l, m, cos(θ), ϕ)
            ψ²[counter] = abs(ψ)^2
            # this is the probability of finding an electron of the orbital (n,l,m) at the position (r,θ,ϕ)
            # transforming spherical coordinates to cartesian coordinates and storing them
            x[counter] = r*sin(θ)*cos(ϕ)
            y[counter] = r*sin(θ)*sin(ϕ)
            z[counter] = r*cos(θ)
            counter = counter +1
        end
    end
end
##
ψ²_sum =sum(ψ²) # This should be close to 1
ψ² = ψ²/ψ²_sum # Now it is normalized
# Initialize the points that need to be plotted
X = []
Y = []
Z = []
c=0
for b=1:PH
    if ψ² < 1.1 && ψ² > 0.9
        c = c + 1
        push!(X,x[c])
        push!(Y,y[c])
        push!(Z,z[c])
            # I am storing all the points at which 0.9<Psi/Mpsi<1.1 
            # Sometimes I changed the range to get better plot
    b=b+1

# Now plotting the stored points
Plots.scatter3d(X, Y, Z, markersize = 2, color = :red)
title!("n= l=2 m=2")
xlabel!("x")
ylabel!("y")
zlabel!("z")

I am new to Julia and I have just migrated this MATLAB code to Julia.
I keep getting the following error when running the above code

┌ Warning: Assignment to `counter` in soft scope is ambiguous because a global variable by the same name exists: `counter` will be treated as a new local. Disambiguate by using `local counter` to suppress this warning or `global counter` to assign to the existing global variable.
└ @ c:\Users\vikra\OneDrive\Desktop\Julia Programmes\Hydrogen Atom\Hydrogen_atom.jl:114
0.0ERROR: StackOverflowError:
Stacktrace:
 [1] _repeat_apply(f::Function, n::Int64) (repeats 43392 times)
   @ Symbolics C:\Users\vikra\.julia\packages\Symbolics\OrNx6\src\diff.jl:366
 [2] _repeat_apply (repeats 2 times)
   @ C:\Users\vikra\.julia\packages\Symbolics\OrNx6\src\diff.jl:366 [inlined]
 [3] ^
   @ C:\Users\vikra\.julia\packages\Symbolics\OrNx6\src\diff.jl:46 [inlined]
 [4] AssoLegrande(l::Int64, m::Int64, cst::Float64)
   @ Main c:\Users\vikra\OneDrive\Desktop\Julia Programmes\Hydrogen Atom\Hydrogen_atom.jl:19
 [5] Angular(l::Int64, m::Int64, cst::Float64, phi::Float64)
   @ Main c:\Users\vikra\OneDrive\Desktop\Julia Programmes\Hydrogen Atom\Hydrogen_atom.jl:75
 [6] top-level scope
   @ c:\Users\vikra\OneDrive\Desktop\Julia Programmes\Hydrogen Atom\Hydrogen_atom.jl:107

Please Help me get the code running

EDIT: sorry, this is only addressing the soft scope warning. It’s worth fixing that, but it isn’t your actual error.

The error you are seeing is telling you that there is a counter outside the loop and one inside the loop and it isn’t certain whether the inside one should be its own thing or the same as the outside.

“Soft scope” is a somewhat messy compromise that helps Julia code run how you would expect on the command line (outside of functions). You can google it to find more discussion, but the short of it is that code at “top level” (not within a function or other construct) in the REPL has slightly unusual scoping rules and tends to create these sorts of issues.

Note that the global in your global counter = 1 initialization does nothing, because all variables outside of other scopes (mostly functions) are global in their module (the REPL is also a module and its name is Main). Note that you didn’t need to label l, m, or n as global. The reason they didn’t have trouble because they weren’t assigned within the loop, so they must have come from outside. One might contend that counter was read before it was written within the loop, so also must be external, but I assume it’s difficult for the compiler to recognize that.

If you add global counter at the top of your ϕ loop, you will tell it that “I want counter to refer to the outer counter and not to a local loop variable.”

Alternatively, this should also be fixed if you put that part of the code within a function, but that’s probably a bigger adjustment.

1 Like

As for your actual error, it would appear it can be reduced to this MWE:

using Symbolics
@variables v
Differential(v)^0

It looks like the 0th differential isn’t supported by Symbolics? If it isn’t but it should be, consider checking whether there is an open issue on its github and opening one if not.

Could you replace this with D_cst_m = iszero(m) ? Returns(cosθ) : Differential(cosθ)^m (and similar for the l version in another function)? Or is there another way you can special case the Differential()^0 case to work?

Anyway, that appears to address the error at hand.

2 Likes

Thank you for your reply. I have wrapped the for loops in a function as that seemed more logical.

I think this was the problem I was running into previously. I have added some if statements to mitigate this

It looks like the 0th differential isn’t supported by Symbolics?

Now fixed Add support for Differential(...)^0 by DanielVandH · Pull Request #1164 · JuliaSymbolics/Symbolics.jl · GitHub

2 Likes