restart; with(plots):

Newton-Raphson

g := x-> 3*tanh(x)-x^3;

plot(g(x), x = -3 .. 3);

x-g(x)/(diff(g(x), x));

f := unapply(%, x);

x1 := 1.2; N:=20;

xx := Array(1..N):

xx[1] := x1:

for n from 1 to N-1 do

xx[n+1] := evalf( f(xx[n]) );

end do:

pp := seq([i, xx[i]], i = 1 .. N):

pointplot([pp], labels = (["n", "xn"]));

xstar := pp[N][2];

check with maple's algorithm

fsolve(g(x), x = 1);

Not bad..

df := D(f);

df(xstar);

the stability is 0, i.e. superstable; this means that the series converges extremely fast

g := 'g';

f := x-g(x)/(diff(g(x), x));

diff(f, x);

Since g(xstar) = 0, we conclude (naively) that df/dx(xstar) = 0 for arbitrary g

i.e. the system is always superstable whatever g.

....which is incorrect, as will follow below

g := x-> (cosh(4-3/x)-1)^(1/a);

x-g(x)/(diff(g(x), x));

f := unapply(%, x);

df := D(f);

a := 6;

plot(g(x), x = .7 .. .8, labels = (["x", "g(x)"]));

plot([x, f(x)], x = 0 .. 1, labels = (["x", "f(x)"]), legend = (["y=x", "y=f(x)"]));

x1 := 0.4; N:=200;

xx := Array(1..N):

xx[1] := x1:

for n from 1 to N-1 do

xx[n+1] := evalf( f(xx[n]) );

end do:

pp := seq([i, xx[i]], i = 1 .. N):

pointplot([pp], labels = (["n", "xn"]));

ret := seq([xx[i], xx[i+1]], i = 1 .. N-1):

pointplot([ret], labels = (["xn", "xn+1"]));

apparently a chaotic series is generated

g(3/4);

So there is no question that xstar = 3/4 is the right solution. However:

f(3/4);

df(3/4);

Lets try the limit

limit(f(x), x = 3/4);

limit(df(x), x = 3/4);

We see that df/dx /= 0, as concluded before.

Bifurcation diagram

a := 'a';

N := 200; Nr := 100; x1 := 0.1; amin := 3; amax := 6;

Nmin := N/2; Nmax := N;

Define Arrays to store the bifurcation points and r-values

Npts := Nmax-Nmin+1;

points := Array(1..Nr, 1..Npts);

ri := Array(1..Nr);

Calculate x[n] for various r

xx := Array(1..N): xx[1] := x1:

for ir from 1 to Nr do

a := evalf(amin + (ir-1) / (Nr-1) * (amax - amin));

ri[ir] := a;

Calculate the x[n] series

for n from 1 to N-1 do

xx[n+1] := evalf( f(xx[n]) );

end do:

Add the converged last section of the series to the bifurcation list

points[ir, 1..Npts] := xx[Nmin..Nmax]:

end do:

Plot the bifurcation diagram. Firstly, we 'serialize' the points Array, by creating an pp Array which contains [r, x] pairs.

pp := Array(1..Nr*Npts, 1..2):

for ir from 1 to Nr do

pp[1 + (ir-1)*Npts..ir*Npts, 1] := ri[ir];

pp[1 + (ir-1)*Npts..ir*Npts, 2] := points[ir];

end do:

Convert the Array to a list and plot the bifurcation diagram

convert(pp, listlist):

pointplot(%, symbol=point, symbolsize=1,labels=["a", "x"]);

For a = 4.5 there seems to be a period 2

a := 4.5;

f2 := x-> f(f(x));

f2(x):

plot([f(x), f2(x), x], x = .2 .. 1);

x2a := fsolve(f2(x) = x, x = .6);

x2b := f(x2a);

df2 := D(f2):

df2(x2a);

df2(x2b);

sigma := df(x2a)*df(x2b);

Three ways that show that x2a,x2b form a stable period -2 solution, since the stability |sigma| < 1

Will maple be fooled? Lets try:

a := 6;

fsolve(g(x), x = .4);

Clearly not.