Thanks for everybody's help. I got something working which satisfies
me.
Regarding Cedric's question, I am trying to compute the Feigenbaum
number delta to arbitrary high precision.
https://en.wikipedia.org/wiki/Feigenbaum_constants
I am using a brute-force method to compute delta. The basic idea is
that each bifurcation point may be found by solving a polynomial
system of order N+1, where N is the number of periods. Delta is then
found in the limit via a ratio of bifurcation points. I'm using
Newton's method to solve the system giving the bifurcation points. My
method is not particularly clever, but it's direct and easily
understandable for a math-dud like me.
Julia has two features which make this computational method possible:
1. I can generate code on the fly using metaprogramming. That's how I
generate a function implementing my N+1st order system (and the
associated Jacobean) which I can then call from Newton's method
at each bifurcation point. 2. Ability to easily use BigFloat as my
variable type, giving aribtrary precision. I do a lot of work in
Matlab, and you can't do either of these things in Matlab.
Regarding generating the code, I fiddled around with directly
generating the AST briefly. Then I decided that it was a little too
tricky for my knowledge of Julia at this point. I decided that I
liked the approach where I generate a string, and then parse and eval
it. The reason is that I can easily inspect and verify a textual
representation of my program since it is simply human-readable Julia
code. The AST is one level of indirection away, and I'm not an expert
enough to be able to read it as easily as Julia code itself.
Here's what I eventually did. This builds the textual representation
of my polynomial system:
function build_f(N)
mymod(x, y) = mod(x-1, y)+1
P = string("function f(X::Vector)\n")
P = string(P, "z = zeros(X);\n")
for i in 1:N
P = string(P, "x$i = X[$i];\n");
end
P = string(P, "mu = X[$(N+1)];\n")
for i in 1:N
ip1 = mymod(i+1, N)
P = string(P, "z[$i] = mu*x$i*(1-x$i) - x$ip1;\n")
end
Q = "mu"
for i = 2:N
Q = string(Q, "*mu")
end
P = string(P, "mu$N = $Q;\n")
for i in 1:N
P = string(P, "xx$i = 1-2*x$i;\n");
end
Q = "mu$N"
for i in 1:N
Q = string(Q, "*xx$i")
end
P = string(P, "t = $Q;\n")
P = string(P, "z[$(N+1)] = t*t-1;\n")
P = string(P, "return z;\n");
P = string(P, "end;\n")
return P
end
And this is the result of calling it for N=4:
julia> y = build_f(4)
"function f(X::Vector)\nz = zeros(X);\nx1 = X[1];\nx2 = X[2];\nx3 =
X[3];\nx4 = X[4];\nmu = X[5];\nz[1] = mu*x1*(1-x1) - x2;\nz[2] =
mu*x2*(1-x2) - x3;\nz[3] = mu*x3*(1-x3) - x4;\nz[4] = mu*x4*(1-x4) -
x1;\nmu4 = mu*mu*mu*mu;\nxx1 = 1-2*x1;\nxx2 = 1-2*x2;\nxx3 =
1-2*x3;\nxx4 = 1-2*x4;\nt = mu4*xx1*xx2*xx3*xx4;\nz[5] =
t*t-1;\nreturn z;\nend;\n"
julia> println(y)
function f(X::Vector)
z = zeros(X);
x1 = X[1];
x2 = X[2];
x3 = X[3];
x4 = X[4];
mu = X[5];
z[1] = mu*x1*(1-x1) - x2;
z[2] = mu*x2*(1-x2) - x3;
z[3] = mu*x3*(1-x3) - x4;
z[4] = mu*x4*(1-x4) - x1;
mu4 = mu*mu*mu*mu;
xx1 = 1-2*x1;
xx2 = 1-2*x2;
xx3 = 1-2*x3;
xx4 = 1-2*x4;
t = mu4*xx1*xx2*xx3*xx4;
z[5] = t*t-1;
return z;
end;
I'm very pleased I can inspect the code in this format and make sure I
haven't made any syntax errors.
To solve the system, I pass the fcns to Newton's method like this:
# Create fcn whose roots we want to find.
y = build_f(N);
Pf = eval(parse(y));
# Build Jacobean of fcn.
y = build_J(N);
PJ = eval(parse(y));
# Do N dimensional Newton's method.
Xstar = newtonND(f, J, X0, tol);
I probably don't need to assign the result of eval(parse()) to a
variable Pf since eval(parse()) drops the name of the function into my
current namespace, right?
The other question is: does this approach create performance-degrading
global objects? If so, how do I mitigate that problem?
At this point I am struggling with getting Newton's method to converge
properly for high N. However, that has to do with selecting the
correct initial condition -- i.e. a typical numerical analysis problem
-- and is unrelated to the specifics of Julia. What's cool is that I
have spent only 1 1/2 days on this project, and am quite far along on
completing it. If I had to stitch everything together using, say, C++
or Mathematica, I'd probably still be struggling with the question of
how to create the system I need to solve.
Stuart
On Sat, 26 Dec 2015, Cedric St-Jean wrote:
That is, I'm solving a set of N polynomial equations, where N is a
parameter.
That sounds interesting, could you post some more details? At this stage,
it's not clear if you really need to generate a new program at all. But if
you do want to do that, then for sure, building the AST will be easier and
more reliable than generating text.
Best,
C?dric
On Saturday, December 26, 2015 at 3:24:51 PM UTC-5, Stuart Brorson wrote:
Cool. Thanks for the code example. I will play with it to see what I
can do with it.
Thanks,
Stuart
On Sat, 26 Dec 2015, Josh Langsfeld wrote:
As Yichao is hinting, you may find a macro to be a cleaner way of making
your functions instead of constructing and parsing a string.
macro return_fcn(N)
xexprs = Expr[:($(symbol(:x,i)) = X[$i]) for i=1:N]
return esc(:(
function $(symbol(:f,N))(X::Vector)
$(xexprs...)
mu = X[$(N+1)]
end
))
end
julia> macroexpand(:(@return_fcn(2)))
:(function f2(X::Vector)
x1 = X[1]
x2 = X[2]
mu = X[3]
end)
julia> macroexpand(:(@return_fcn(4)))
:(function f4(X::Vector)
x1 = X[1]
x2 = X[2]
x3 = X[3]
x4 = X[4]
mu = X[5]
end)
On Saturday, December 26, 2015 at 2:03:56 PM UTC-5, Stuart Brorson
wrote:
Julia users,
I'm fiddling around with Julia's strings & metaprogramming. I am
constructing a function by concatenating a bunch of strings together
to create my function, like this:
function return_fcn(N)
P = string("function f$N(X::Vector)\n")
for i in 1:N
P = string(P, "x$i = X[$i];\n");
end
P = string(P, "mu = X[$(N+1)];\n")
etc....
When I execute this code, I get:
julia> y = return_fcn(2)
"function f2(X::Vector)\nx1 = X[1];\nx2 = X[2];\nmu = X[3];\n"
However, what I really want to see is
function f2(X::Vector)
x1 = X[1];
x2 = X[2];
mu = X[3];
"show(y)" doesn't seem to do what I want. Later, when I do
eval(parse(y))
then I get a function which executes correctly. My problem is simply
that I can't get Julia to give me a string I can read easily. This
will be a very big issue for me when N -> 1024, 2048, etc....
Questions:
1. How can I escape the \n to get a real <CR><LF> in my displayed
string?
2. Is this the optimal way to construct a program for later execution
(i.e. metaprogramming)?
Thanks for all wisdom you have to offer.
Stuart
