# Lesson 8 # ====== # # Problem: find a cubic polynomial for which Newton's method has an attracting 2-cycle. > restart; > newt:= x -> x - f(x)/D(f)(x); > dnewt:=D(newt); # Since the derivative of newt is 0 when f''(x)=0, the simplest way to ensure that a cycle is # attracting is to make this be true at one point on # the cycle. # Here is a cubic polynomial. > f:= x -> x^3 + a*x^2 + b*x + c; # We choose the points on the cycle to be 0 and 1. The first two equations make this a # cycle. > eq1:= newt(0) = 1; > eq2:= newt(1) = 0; # The third equation is f''(0)=0. > eq3:= (D@@2)(f)(0) = 0; # "solve" will solve (if possible) a set of equations for a set of variables. > s:=solve({eq1, eq2, eq3}, {a,b,c}); # Here's our polynomial. > subs(s, f(x)); # We can use "subs" on a function definition as well as on an expression. # However, we must use "eval" on the name of the function, otherwise we'd just get the # name f and not the definition. > f:= subs(s, eval(f)); # Now to try a picture of this. The "newtstep" command plots one step # of Newton's method. > newtstep:= x -> plot([ [x,0], [x, f(x)], [newt(x), 0], [newt(x), f(newt(x))] ]); # We choose a starting point near, but not on, the cycle. Some experimenting was # necessary: 0.2 would not have worked. > x[0]:= 0.1; # We calculate 10 steps of Newton's method. > for jj from 1 to 10 do x[jj]:= newt(x[jj-1]); od: > with(plots): # We plot these steps, together with the graph of f. > display({ seq( newtstep(x[jj]), jj=0 .. 10), > plot(f(x), x = -2 .. 1.5) }); # Here are the even-numbered x[j]'s. Note that they approach 0, at # first slowly but then very rapidly. The quadratic convergence to 0 # (x[j+2] approximately k x[j]^2) is a result of the fact that D(newt) is # 0 at one of the points of the cycle. > seq(x[2*jj], jj=0..5); # How big is the "basin of attraction" of this 2-cycle? > plot({-1,1,newt(newt(x))/x}, x = -1 .. 1, -2 .. 2); # If x[0] is between about -.12 and .14, x[2] will be closer to 0 than x[0] # was, and the sequence x[j] will approach the 2-cycle. Outside this # interval, we can't be sure what will happen. # # # "solve" in Several Variables # ==================== > restart; # Here are two polynomials in two variables. > p1:= (x,y) -> x^4 - x*y^2 - y^3 + 1 ; > p2:= (x,y) -> x^3 - 5* x^2 * y + y^4 - 2 ; # Do you think Maple can find a solution for the simultaneous equations p1(x,y) = 0 and # p2(x,y)=0? > sol:=solve({p1(x,y)=0, p2(x,y)=0}, {x,y}); # It did find the solutions after a fashion: x is "RootOf" a complicated polynomial, and y is a # complicated polynomial of that. It's our bad luck that the first complicated polynomial # can't be solved in "closed form". But what is this polynomial? > res:=resultant(p1(x,y),p2(x,y),y); # This is the "resultant" of p1(x,y) and p2(x,y), considered as polynomials in y (with # coefficients that are polynomials in x). The main property is that if two polynomials in y # have a common root, their resultant is 0. # # To understand resultants, it's simplest to eliminate the distraction of the x variable, and # start with polynomials with constant coefficients. I'll take these same polynomials, but fix # x=3. > f:= unapply(p1(3,y), y); > g:= unapply(p2(3,y), y); # Now the resultant of f(y) and g(y) is a number. > resultant(f(y),g(y),y); # Here's the definition of resultant. Suppose f(x) = a x^n + ... and # g(x) = b x^m + ... are polynomials of degrees n and m respectively. # Let f(x) have roots r[1], ..., r[n] (with perhaps multiple roots) and g(x) have roots s[1], ..., # s[m]. Then the resultant of f and g is > a^m * b^n * Product(Product(r[i] - s[j], i = 1 ... n), j = 1 ... m); # Note that this will be 0 if f and g have a root in common. Let's calculate the roots for our # examples. > r:= fsolve(f(y), y, complex); > s:= fsolve(g(y), y, complex); > product(product(r[i]-s[j], i=1..3), j=1..4); > round("); # Here's a shortcut to computing it. Note that > 'g'(x) = b * Product(x-s[j], j = 1 .. m); # So the resultant of f and g must be > a^m * Product('g'(r[i]), i = 1 .. n); # Let's try it. > product(g(r[i]), i=1..3); # It also works the other way, write > 'f'(x) = a * Product(x - r[i], i = 1 .. n); # but note that each factor (r[i] - s[j]) must be multiplied by -1, so the resultant is > (-1)^(m*n) * b^n * Product('f'(s[j]), j = 1 .. m); > product(f(s[j]), j=1..4); # Actually, there is a special form of "product" for a product of some expression over the # roots of a polynomial. It doesn't resort to floating-point approximations: it turns out # there are efficient ways to calculate the sum of the k-th powers of the roots of a # polynomial, where k is any positive integer. For example, the sum of the roots is # - the coefficient of x^(d-1), where d is the degree. > product(f(t), t = RootOf(g(y)) ); # Let's try this for our original polynomials, which are functions of x and y, thinking of # them as polynomials in y whose coefficients happen to involve x. > product(p1(x,t), t = RootOf(p2(x,y), y) ); > " - res; # If you had to compute the roots of one of the polynomials to calculate the resultant, it # wouldn't be very useful. Fortunately, there is an efficient way to calculate the resultant. # It is based on the principle that for any polynomial p, the resultant of f and g + p f is (up # to a power of a) the resultant of f and g. That's because # (g + p f)(r[i]) = g(r[i]). # Here's how we could use this to calculate the resultant of our f and g by hand. # We'd start with "long division" of g by f. > 'g(y)' = quo(g(y),f(y), y) * 'f(y)' + rem(g(y), f(y), y); > normal("); # So resultant(f(y), g(y)) = resultant(f(y), -221 + 37 y + 9 y^2). This reduces the problem # from calculating the resultant of a quartic and a cubic to a cubic and a quadratic. # Repeating the process, we get a quadratic and a linear polynomial, and then linear and # constant. > g1:= rem(g(y), f(y), y); > resultant(f(y), g(y),y) = resultant(f(y), g1,y); > f1:= rem(f(y), g1, y); > resultant(f1, g1, y) = resultant(f(y), g1, y)/9^2; > g2:= rem(g1,f1,y); > resultant(f1,g2,y)= resultant(f1,g1,y)*(-2359/81)^(-2); >