## 192 Reputation

15 years, 256 days

Thanks for that!

## Thanks for that!...

Thanks for that!

g is fine but the last line should be...

seq(a[2*i]=coeff(g(x,2*N + 2),x^(2*i)),i=1..N);

Otherwise all the coefficients equate to zero.

Also, the line,

solve({%},{seq(a[2*i],i=1..N)},AllSolutions=true);

Appears to work in that I can see what it's doing but it's not giving the correct solutions.

Would fsolve not be a better option for this type of question? Even when N=3 solve takes 47 seconds to run using 55Mb of memory while my fsolve one runs in around 0.14 seconds using hardly any memory.

I might have to give up on this one since having spoken to a few of the lecturers around here who work with Maple, they don't think this will work for any high degree polynomial. This is basically a very inneficient way of finding feigenbaums constants but I wanted to see what I could do.

g is fine but the last line should be...

seq(a[2*i]=coeff(g(x,2*N + 2),x^(2*i)),i=1..N);

Otherwise all the coefficients equate to zero.

Also, the line,

solve({%},{seq(a[2*i],i=1..N)},AllSolutions=true);

Appears to work in that I can see what it's doing but it's not giving the correct solutions.

Would fsolve not be a better option for this type of question? Even when N=3 solve takes 47 seconds to run using 55Mb of memory while my fsolve one runs in around 0.14 seconds using hardly any memory.

I might have to give up on this one since having spoken to a few of the lecturers around here who work with Maple, they don't think this will work for any high degree polynomial. This is basically a very inneficient way of finding feigenbaums constants but I wanted to see what I could do.

## Cool thanks for that!   I had a lit...

Cool thanks for that!

I had a little play around with the code... Couldn't get the answers I was after but some things I noted were that...

f contains x terms and g doesn't so f - g doesn't really work.   yet...

alpha should probably be alpha := n-> (f(1,n/2))^(-1): or even alpha := (f(1,N))^(-1):

Similarly g := (x,n) -> simplify(coeff(series(alpha(n)*f(f(x/alpha(n),n),n),x), x^n)): should be

g := (x,n) -> simplify(coeff(series(alpha*f(f(x/alpha,N),N),x), x^n)):

(I think...! and also using my second definition of alpha.)

I know you've not had time to look at it so I realize it's not the finished thing but that's just what I noticed.

Thanks for the code so far though. I'll look more at it however I'm beginning to think that this will just involve equations that are too large no matter how efficient my code is...

## Cool thanks for that!   I had a lit...

Cool thanks for that!

I had a little play around with the code... Couldn't get the answers I was after but some things I noted were that...

f contains x terms and g doesn't so f - g doesn't really work.   yet...

alpha should probably be alpha := n-> (f(1,n/2))^(-1): or even alpha := (f(1,N))^(-1):

Similarly g := (x,n) -> simplify(coeff(series(alpha(n)*f(f(x/alpha(n),n),n),x), x^n)): should be

g := (x,n) -> simplify(coeff(series(alpha*f(f(x/alpha,N),N),x), x^n)):

(I think...! and also using my second definition of alpha.)

I know you've not had time to look at it so I realize it's not the finished thing but that's just what I noticed.

Thanks for the code so far though. I'll look more at it however I'm beginning to think that this will just involve equations that are too large no matter how efficient my code is...

## Ok no worries thanks anyway.   I'...

Ok no worries thanks anyway.

I'll give it one more explanation and if still no replies I guess I'm doomed to solve this alone...

Ignore everything I've said and just concentrate on this.

How do you solve this equation?

(a + b)x + (a - b) = 3x + 4?

You compare coefficients!

We have,

(a + b) = 3 and

(a - b) = 4

This should give a = 3.5 and b = 0.5.

this is basically what I'm trying to do. Compare coefficients from two functions to find the value of the constants.

In my case, the functions are,

f(x) = 1 - lambda*x^2

and

alpha*f(f(x/alpha))

where alpha = 1/f(1) = 1/(1 - lambda)

Now the thing is, when you substitute f(x) = 1 - lambda*x^2 into alpha*f(f(x/alpha)) you get a function which also has x^4 terms which we do not have in f(x). So we just discard and ignore these higher order terms.

So working with the functions as above, we get,

alpha*f(f(x/alpha)) = 1 + (2*lambda^2 - 2*lambda^3)*x^2 + terms of order x^4 which we ignore.

So now we compare this to f(x) = 1 - lambda*x^2.

We then see that in order for the two functions to be equal (again, ignoring the x^4 terms), we must solve,

(2*lambda^2 - 2*lambda^3) = - lambda

This gives us lambda = (1 + sqrt(3))/2

We then substiture this back into out definition of alpha to get alpha being roughly -2.7.

This is actually a known constant (one of Feigenbaums constants) for which the true value is -2.5... so we are close but can improve.

Now! We add on extra terms to f(x) and repeat the above proceedure except now we will have more equations to solve with more unknowns.

So if we let f(x) = 1 - lambda*x^2 + a*x^4.

We now have alpha = 1/f(1) = 1/(1 - lambda + a)

Now if we were to write out alpha*f(f(x/alpha) in full, we'd see we have terms up x^16. We then discard anything above our highest term in f(x) which is x^4.

Now we will have coefficients in alpha*f(f(x/alpha)) we can compare to coefficients in f(x) and solve for a and lambda and thus find a value for alpha.

So! That's what I'm doing.

The more terms you add onto f(x), the more accurate alpha will get.

This is the code I'm currently using

with(RealDomain):
Order := 12:
f := x -> 1 - lambda*x^2 + a*x^4 + b*x^6 + c*x^8 + e*x^10:
alpha := 1/(1 - lambda + a + b + c + e):
Lambda := simplify(coeff(series(alpha*f(f(x/alpha)),x),x^2)):
A := simplify(coeff(series(alpha*f(f(x/alpha)),x), x^4)):
B := simplify(coeff(series(alpha*f(f(x/alpha)),x), x^6)):
C := simplify(coeff(series(alpha*f(f(x/alpha)),x), x^8)):
E := simplify(coeff(series(alpha*f(f(x/alpha)),x), x^10)):
Eq := {Lambda = -lambda,A = a,B = b,C = c,E = e}:
fsolve(Eq,{lambda,a,b,c,e},avoid={lambda=0,a=0,b=0,c=0,e=0});

But it will only go as high as x*10. Have left on a version where I had x^12 as the highest term but never computed even after a day of running.

So, if you understood this... Any idea how I can make my code more efficient.

## Sorry, another point to make.   T...

Sorry, another point to make.

The above code, while it seems to work... Isn't giving me the roots I was expecting however, it does give me the correct roots here when I'm just having one extra a*x^4 term. The third solution is the one I'm after.

with(RealDomain):
Order := 6:
f := x -> 1 - lambda*x^2 + a*x^4:
alpha := 1/(1 - lambda + a):
Lambda := simplify(coeff(series(alpha*f(f(x/alpha)),x),x^2)):
A := simplify(coeff(series(alpha*f(f(x/alpha)),x), x^4)):
with(SolveTools):
PolynomialSystem([Lambda + lambda,A - a], [lambda,a], 100):
evalf(%);

## I've had 2 emails saying there has been ...

I've had 2 emails saying there has been posts here but I can't see them so...

But anyway, here's my updated code. Perhaps this will make it clearer what I'm trying to do.

with(RealDomain):
Order := 8:
f := x -> 1 - lambda*x^2 + a*x^4 + b*x^6:
alpha := 1/(1 - lambda + a + b):
Lambda := simplify(coeff(series(alpha*f(f(x/alpha)),x),x^2)):
A := simplify(coeff(series(alpha*f(f(x/alpha)),x), x^4)):
B := simplify(coeff(series(alpha*f(f(x/alpha)),x), x^6)):
with(SolveTools):
PolynomialSystem([Lambda + lambda,A - a,B - b], [lambda,a,b], 100):
evalf(%);

Note I had 100 in PolynomialSolve since I wasn't sure how many solutions I would come across.

Most of the solutions found have one of the terms being zero. Is there anyway I can stop this happening? I.e. all terms are non zero.

## Righto thanks for the reply.   I ...

I guess I must have mistyped something when I typed in up here...

This is just copied straight from the worksheet. Final answer should be alpha := -2.534...

f := x -> 1 - lambda*x^2 + a*x^4:
g := 1 - lambda*x^2 + a*x^4:
alpha := 1/f(1):
Lambda := select(has, collect(alpha*f(f(x/alpha)),x), x^2):
lambda := [solve(Lambda = select(has, g, x^2), lambda)][2];

A := select(has, collect(alpha*f(f(x/alpha)),x), x^4):
a := [evalf(solve(A = select(has, g, x^4), a))][1];
alpha := evalf(1/f(1));

The idea behind the code is to take two equations, one being comparing the coefficients from the x^2 terms and finding lambda as a function of a.

Then compare the coefficients of the x^4 terms which should give an expression involving just a. This can then be solved and hence lambda then alpha can be found.

## Righto thanks for the reply.   I ...

I guess I must have mistyped something when I typed in up here...

This is just copied straight from the worksheet. Final answer should be alpha := -2.534...

f := x -> 1 - lambda*x^2 + a*x^4:
g := 1 - lambda*x^2 + a*x^4:
alpha := 1/f(1):
Lambda := select(has, collect(alpha*f(f(x/alpha)),x), x^2):
lambda := [solve(Lambda = select(has, g, x^2), lambda)][2];

A := select(has, collect(alpha*f(f(x/alpha)),x), x^4):
a := [evalf(solve(A = select(has, g, x^4), a))][1];
alpha := evalf(1/f(1));

The idea behind the code is to take two equations, one being comparing the coefficients from the x^2 terms and finding lambda as a function of a.

Then compare the coefficients of the x^4 terms which should give an expression involving just a. This can then be solved and hence lambda then alpha can be found.

## Yes g is just f but f has the x -> pa...

Yes g is just f but f has the x -> part which g doesn't have. I couldn't get the select(has, ...) part to work with f since I had the x -> ... part.

lambda doesn't depend on A. The Lambda and lambda part writes lambda as a function of a (and b for the second part).

Then the A := ... and a := then solve for what a is after subing in the lambda which will be a function of a.

I'll have to go through the steps I see but I thought it was fairly simple to figure out what was trying to be achieved.

This is all based on the Feigenbaum-Cvitanovi´c equation for period-doubling and about finding a value for the universal constant alpha.

Take the function f(x) := 1 - lambda*x^2

alpha is defined (in this case) to be 1/f(1) which is 1/(1-lambda)

Now we use Feigenbaum-Cvitanovi´c equation which is;

alpha(f(f(x/alpha)) = f(x) to find a value for lambda and hence alpha.

However, simply subing f(x) into this doesn't give an exact solution, i.e. for f(x) = 1 - lambda*x^2, alpha(f(f(x/alpha)) does not exactly equal f(x).

However we can find a value of lambda such that...

alpha(f(f(x/alpha)) = f(x) + O(x^4)

where O is the big O notation for order.

The first line is f(x)

the last line is alpha*f(f(x/alpha)).

By comparing coefficients and treating the x^4 part as a higher order term and ignoring it, we see that we must find lambda such that,

2*lambda^2*(-lambda + 1) = -lambda

i.e. the coefficients of x^2 from each equation.

Doing this gets a lambda of (1 + sqrt(3))/2 which we sub into the the equation for alpha (1/f(1)) which gives alpha ~ -2.7

With me so far?

Then now I want to improve on this result.

This is done by adding more terms to f(x).

So we set f(x) = 1 - lambda*x^2 + a*x^4.

Now if we sub this into the Feigenbaum-Cvitanovi´c equation we will have two unknowns, lambda and a, and two equations with which to find them.

I.e......

We will get alpha*f(f(x/alpha)) = 1 + x^2*(something) + x^4*(something) + O(x^6)

With which we can compare it to f(x), disregard the O(x^6) term and compare coefficients to find values for lambda and a then finally, alpha.

This results in alpha = -2.53... which is much closer to the actual value of alpha (-2.50...) than the -2.7 I got earlier.

So now I want to add a second added term, namely +b*x^6 and perform all the above steps again to hopefully get a more accurate answer for alpha.

But I can't get it to work as the equations are just to complex or my code is just not very effiecint.

If you have any specific points you are not understanding could you point them out and I can explain them in more detail. I'd rather avoid getting deep into the maths... My code I posted in the first post DOES work as I have used it with different example but I'm not very good at writing effiecient code.

## Yes g is just f but f has the x -> pa...

Yes g is just f but f has the x -> part which g doesn't have. I couldn't get the select(has, ...) part to work with f since I had the x -> ... part.

lambda doesn't depend on A. The Lambda and lambda part writes lambda as a function of a (and b for the second part).

Then the A := ... and a := then solve for what a is after subing in the lambda which will be a function of a.

I'll have to go through the steps I see but I thought it was fairly simple to figure out what was trying to be achieved.

This is all based on the Feigenbaum-Cvitanovi´c equation for period-doubling and about finding a value for the universal constant alpha.

Take the function f(x) := 1 - lambda*x^2

alpha is defined (in this case) to be 1/f(1) which is 1/(1-lambda)

Now we use Feigenbaum-Cvitanovi´c equation which is;

alpha(f(f(x/alpha)) = f(x) to find a value for lambda and hence alpha.

However, simply subing f(x) into this doesn't give an exact solution, i.e. for f(x) = 1 - lambda*x^2, alpha(f(f(x/alpha)) does not exactly equal f(x).

However we can find a value of lambda such that...

alpha(f(f(x/alpha)) = f(x) + O(x^4)

where O is the big O notation for order.

The first line is f(x)

the last line is alpha*f(f(x/alpha)).

By comparing coefficients and treating the x^4 part as a higher order term and ignoring it, we see that we must find lambda such that,

2*lambda^2*(-lambda + 1) = -lambda

i.e. the coefficients of x^2 from each equation.

Doing this gets a lambda of (1 + sqrt(3))/2 which we sub into the the equation for alpha (1/f(1)) which gives alpha ~ -2.7

With me so far?

Then now I want to improve on this result.

This is done by adding more terms to f(x).

So we set f(x) = 1 - lambda*x^2 + a*x^4.

Now if we sub this into the Feigenbaum-Cvitanovi´c equation we will have two unknowns, lambda and a, and two equations with which to find them.

I.e......

We will get alpha*f(f(x/alpha)) = 1 + x^2*(something) + x^4*(something) + O(x^6)

With which we can compare it to f(x), disregard the O(x^6) term and compare coefficients to find values for lambda and a then finally, alpha.

This results in alpha = -2.53... which is much closer to the actual value of alpha (-2.50...) than the -2.7 I got earlier.

So now I want to add a second added term, namely +b*x^6 and perform all the above steps again to hopefully get a more accurate answer for alpha.

But I can't get it to work as the equations are just to complex or my code is just not very effiecint.

If you have any specific points you are not understanding could you point them out and I can explain them in more detail. I'd rather avoid getting deep into the maths... My code I posted in the first post DOES work as I have used it with different example but I'm not very good at writing effiecient code.

 1 2 Page 2 of 2
﻿