Hi,
I'm trying to solve system of nonlinear equations ( e1 and e2) to find c
value.
Below is a program which solves it for period 1, 2 and some 3 ( see
Curve3a).
For other values gives error :
Curve3b:GiveListAngle(3,1,center3[2],50)
"mnewton: the process doesn't converge or it converges too slowly."
LIST-REF: invalid subscript: 1
I have tried to increase :
newtonmaxiter:100;
fpprec:100;
newtonepsilon:bfloat(10^(-fpprec/2));
or use mnewton.mac ver 1.1.
Still I cant get a result.
WHat can I change do make it work ?
Regards
Adam
=== info ===========
build_info()$
Maxima version: 5.22.1
Maxima build date: 22:10 11/8/2010
Host type: x86_64-unknown-linux-gnu
Lisp implementation type: GNU Common Lisp (GCL)
Lisp implementation version: GCL 2.6.7
=============== code start ==============
kill(all);
remvalue(all);
/* functions */
/* complex quadratic polynomial */
f(z,c):=z*z+c;
/* iterated function */
fn(p, z, c) :=
if p=0 then z
elseif p=1 then f(z,c)
else expand(f(fn(p-1, z, c),c));
/* Standard polynomial , which roots are periodic z-points of period p
and its divisors */
F(p, z, c) := fn(p, z, c) - z ;
/* Multiplier of periodic z-point */
m(p):=expand(diff(fn(p,z,c),z,1));
GiveC(period,radius, angle, center):=
block(
[e1,e2,s,w, a],
/* first equation defining periodic points */
e1:F(period,z,c),
/* point of unit circle */
a:angle*2*%pi, /* from turns to radians */
w:radius*cos(a) + radius*sin(a)*%i,
/* second equation defining multiplier of periodic point */
e2:m(period)- w,
/* solutio using newton method */
s: mnewton([e1,e2],[c,z],[center,0 ]),
return(rhs(s[1][1]))
)$
GiveListAngle(period,radius, center, NumberOfPoints ):=
block(
[angleStep,angle, ListOfPoints, point],
angleStep:1.0/NumberOfPoints,
ListOfPoints:[],
for n:1 thru NumberOfPoints step 1 do
( angle:n*angleStep,
point:GiveC(period,radius,angle, center),
ListOfPoints:cons([realpart(point),imagpart(point)],ListOfPoints)
),
return(ListOfPoints)
)$
compile(all);
/* ----------- const ------------------------------ */
NumberOf_Points:100$
center1:0$
center2:-1$
center3:[-1.754877666246692760,
-0.12256116687665-0.74486176661974*%i,
-0.12256116687665+0.74486176661974*%i]$
load("mnewton")$
Curve1:GiveListAngle(1,1, center1, NumberOf_Points)$
Curve2:GiveListAngle(2,1, center2, 50)$
Curve3a:GiveListAngle(3,1, center3[1], 10)$
newtonmaxiter:100;
fpprec:100;
newtonepsilon:bfloat(10^(-fpprec/2));
Curve3b:GiveListAngle(3,1, center3[2], 50)$
Curve3c:GiveListAngle(3,1, center3[3], 50)$
============ end code =====================