and compare it to the result obtained by QuadGK. But I’m getting an error that makes no sense to me. Namely I’m getting the error:
UndefVarError: t not defined
Stacktrace:
[1] top-level scope at /data/GitHub/mine/math/julia-scripts/Simpson's_rule_pendulum_integral.jl:56
[2] top-level scope at In[2]:1
and here’s the kicker, t is defined as equalling 0 at line 54. Here’s my script:
# This script essentially finds the time taken for a
# simple pendulum that starts with zero velocity at the
# positive x-axis to reach the negative x-axis (theta=-pi)
# using Simpson's rule
# t is time and y is theta
using Pkg;
Pkg.add("SpecialFunctions")
Pkg.add("QuadGK")
Pkg.add("PyPlot")
using SpecialFunctions;
using QuadGK;
using PyPlot;
# Acceleration rate due to gravity in metres per second
# squared.
g=9.8
# Length of the pendulum in metres.
l=1.0
# Define our integrand
function f(y)
# abs is used below to prevent complex number errors due to the sqrt
return 1/sqrt(abs(2.0*g/l*sin(y)))
end
# Number of steps
N=1000000000;
# How close to theta=0 we start our integration and how
# close to theta=-pi we end our intgration.
# At theta=0,-pi our function becomes undefined, hence
# why we cannot start at exactly theta=0 or end at exactly
# theta=-pi.
tol=1/(25.0*N);
# Integration interval, [y0, yend]
y0=-tol;
yend=-pi+tol;
# Step size
h=(yend-y0)/N;
#
y=y0;
# Our integration function
function Simpson(h,y,i,N)
if i == 1 || i == N+1
return h/3*f(y)
elseif (i % 2) == 1
return 2*h/3*f(y)
else
return 4*h/3*f(y)
end
end
t = 0.0;
# The actual integration
for i=1:N
t = t + Simpson(h,y,i,N);
y = y + h;
end
# Our QuadGK approximation to the integral
T=abs(quadgk(y -> f(y), 0, -pi)[1]);
# Difference between our QuadGK approximation and our
# approximation using Simpson's method
error = abs(T + t)
(the t = t + Simpson(h,y,i,N); line is line 56). Any ideas why Julia doesn’t seem to understand that I have defined t?
@erlebach is entirely correct, just two points to add:
This behaviour will change in the next Julia version
More importantly, working in global scope is bad for performance as the compiler cannot guarantee the type of values used in the loop. It is highly recommended to wrap your code in functions, which would also get you around this problem. The simplest way to see this is to take your existing code and just wrap it in function main() ... end
I’m afraid, as I said in my last reply removing global t from the loop, even within your simpson() function just returns the original error. Namely, I executed this code in Julia:
# This script essentially finds the time taken for a
# simple pendulum that starts with zero velocity at the
# positive x-axis to reach the negative x-axis (theta=-pi)
# using Simpson's rule
# t is time and y is theta
using Pkg;
Pkg.add("SpecialFunctions")
Pkg.add("QuadGK")
Pkg.add("PyPlot")
using SpecialFunctions;
using QuadGK;
using PyPlot;
# Acceleration rate due to gravity in metres per second
# squared.
g=9.8
# Length of the pendulum in metres.
l=1.0
# Define our integrand
function f(y)
# abs is used below to prevent complex number errors due to the sqrt
return 1/sqrt(abs(2.0*g/l*sin(y)))
end
# Number of steps
N=1000000000;
# How close to theta=0 we start our integration and how
# close to theta=-pi we end our intgration.
# At theta=0,-pi our function becomes undefined, hence
# why we cannot start at exactly theta=0 or end at exactly
# theta=-pi.
tol=1/(25.0*N);
# Integration interval, [y0, yend]
y0=-tol;
yend=-pi+tol;
# Step size
h=(yend-y0)/N;
#
y=y0;
# Our integration function
function Simpson(h,y,i,N)
if i == 1 || i == N+1
return h/3*f(y)
elseif (i % 2) == 1
return 2*h/3*f(y)
else
return 4*h/3*f(y)
end
end
function simpson()
t = 0.0;
y = 0.0;
# The actual integration
for i=1:N
# `global` is not valid within a function
t = t + Simpson(h,y,i,N);
y = y + h;
end
end
simpson()
# Our QuadGK approximation to the integral
T=abs(quadgk(y -> f(y), 0, -pi)[1]);
# Difference between our QuadGK approximation and our
# approximation using Simpson's method
error = abs(T + t)
and I received the error:
UndefVarError: t not defined
Stacktrace:
[1] top-level scope at /data/GitHub/mine/math/julia-scripts/Simpson's_rule_pendulum_integral.jl:70
[2] top-level scope at In[4]:1
Hi,
The error is not happening at the same place : now that you have put your loop in a function, you have to return t from your function so that your abs(T+t) can know t (as t is now a local variable inside the simpson function scope.
Thanks, you’re right. I thought it was a different line, but I think I edited my script such that line 70 became simpson() and that lead to my thinking it was that loop that was still generating the error.
For that code to be legal, the compiler would have to determine that the loop body will run at least once. This is a hard problem in general, even though it would be simple in this particular case.
Put something like a = 0 before the loop, and it will work.
Yes, I realized the trick, and use it. However, this implies that the behavior is different from what I thought.
In Julia, one should never assume that loop variables are accessible from outside the loop context. Thank you.
While I realize this isn’t the “Julia Way”. I’ve found that always initializing my variables with local ensures that the variable is defined AND I know what context it’s in. i.e:
function tst()
local a = 0
for i in 1:10
local b = 2
a = 3
end
print("a= ", a)
end
tst()
This question is very thoroughly answered in the relevant section of the manual. I would recommend using it as your primary source to search for answers of questions about keywords of the language, as they surely will be fully explained in the manual.
function tst()
a = 0
# or this, if you do not need any specific initial value for a
# local a
for i in 1:10
b = 2
a = 3 + b
end
print("a= ", a)
end
tst()
julia> tst()
a= 5
Probably we should add that using local b inside the loop allows one to use the name b for a variable inside the loop which repeats the name of a variable in the outer scope:
julia> function tst()
a = 0 ; b = 1
for i in 1:10
local b = 2
a = 3 + b
end
print("a=",a," b=",b)
end
tst (generic function with 1 method)
julia> tst()
a=5 b=1 # b is not changed in the loop
Although we should never do that, as the code might become a mess.
I think (not completely sure) that the fact that b is local to the scope of the loop allows for possible compiler optimizations which would not be possible otherwise. For instance, in my experience, some codes like this in Julia run faster than the almost equivalent translations of Fortran codes, and I this might be one of the explanations, because in Fortran you do not have that exclusive scope for variables in loops.