# 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
#%%
# 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
# Substitute temp = 2*r/(n*a)
ρ_val = 2*r/(n*a)

# Substitute the numeric value of r into the symbolic expression

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

``````

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`?

2 Likes