> restart; with(plots): > f := x -> x^3 - 7*x^2 + 12*x - 5; > df := x -> 3*x^2 - 14*x + 12; > p1 := plot( f(x), x = 0..5): display(p1); > fsolve( f(x) = 0); > ta := x -> f(a) + df(a)*(x-a); > a := 1.3: t1 := plot(ta(x), x = a-1..a+1): > display({t1, p1}); > a := 1.2: t1 := plot(ta(x), x = a-1..a+1): > display({t1, p1}); > a := 1.1: t1 := plot(ta(x), x = a-1..a+1): > display({t1, p1}); > x0 := 1.3: > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := 1.1: > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > peaks := fsolve( df(x) = 0); # Look for an oscillation point. Is there an initial guess a0 with the # property that Newton's method brings us back to a0 after two steps? > a1 := a0 - f(a0)/df(a0); > a2 := a1 - f(a1)/df(a1); > c := fsolve(a2 = a0, a0, a0 = 1..2); > x0 := c; > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > w := c-peaks[1]; > x0 := c - w/2; > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := c - w/4; > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := c - w/8; > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := c - w/16; > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := c - w/32; > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := c - w/64; > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := c - w/128; > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := c - w/256; > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := c - w/10; > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := c - w/100; > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := c - w/1000; > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := c - w/10000; > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); # **************************************************************** # A simpler example of the same phenomenon, one where the oscillation # points can be computed explicitly # **************************************************************** > f := x -> x^3 - x; df := x -> 3*x^2 - 1; > solve(f(x) = 0); plot( f(x), x = -2..2); > x0 := 1/sqrt(5); > x0 := x0 - f(x0)/df(x0); > x0 := 1/sqrt(5.0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := 1/sqrt(5.0) + 1/10^6; > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := 1/sqrt(5.0) + 1/10^7; > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); > x0 := x0 - f(x0)/df(x0); # ********************************************************************** # ****