# Math 210 Sec. 201 # Solutions to Assignment 1 # # 1(a) > eqn:= x^4 - 11*x^3 + 54*x^2 -108*x + 55 = 0; 4 3 2 eqn := x - 11 x + 54 x - 108 x + 55 = 0 > _EnvExplicit:= true; _EnvExplicit := true > exactsoln:= solve(eqn, x); exactsoln := 1/2 1/2 11/4 + 1/12 3 %2 / 1/3 1/2 1/2 2/3 1/2 1/2 1/3\1/2 |414 %1 %2 + 36 %2 %1 + 48 %2 + 3258 3 %1 | + 1/12 I |-------------------------------------------------------------| | 1/3 1/2 | \ %1 %2 / , 1/2 1/2 11/4 + 1/12 3 %2 / 1/3 1/2 1/2 2/3 1/2 1/2 1/3\1/2 |414 %1 %2 + 36 %2 %1 + 48 %2 + 3258 3 %1 | - 1/12 I |-------------------------------------------------------------| | 1/3 1/2 | \ %1 %2 / , 1/2 1/2 11/4 - 1/12 3 %2 / 1/3 1/2 1/2 2/3 1/2 1/2 1/3\1/2 |- 414 %1 %2 - 36 %2 %1 - 48 %2 + 3258 3 %1 | + 1/12 |---------------------------------------------------------------| | 1/3 1/2 | \ %1 %2 / , 1/2 1/2 11/4 - 1/12 3 %2 / 1/3 1/2 1/2 2/3 1/2 1/2 1/3\1/2 |- 414 %1 %2 - 36 %2 %1 - 48 %2 + 3258 3 %1 | - 1/12 |---------------------------------------------------------------| | 1/3 1/2 | \ %1 %2 / 1/2 %1 := 679/2 + 1/18 37343553 1/3 2/3 - 69 %1 + 12 %1 + 16 %2 := -------------------------- 1/3 %1 # Note that "exactsoln" is a sequence. To form a list, we put brackets [] around it. # # (b) > Digits:= 14; Digits := 14 > evalfsoln:= evalf([exactsoln]); evalfsoln := [3.6431728008923 + 3.3827531371481 I, 3.6431728008923 - 3.3827531371481 I, 2.9624959932669, .7511584049487] # (c) > fsolvesoln:= [fsolve(eqn, x)]; fsolvesoln := [.75115840494860, 2.9624959932669] # (d) The second answer from (c) was the same as one of the answers from (b), but the # first was slightly different. To see which was really more accurate (without more # roundoff complicating matters), it's best to increase Digits some more. > Digits:= 20; Digits := 20 > evalfacc:= subs(x=evalfsoln[4], eqn); -11 evalfacc := -.4568114*10 = 0 > fsolveacc:= subs(x=fsolvesoln[1], eqn); -12 fsolveacc := -.188170*10 = 0 # It appears that the solution from (c) was more accurate (although the solution from (b) # was not too far off). The error in (b) is roundoff error from calculating the rather # complicated exact formula for the roots. The reason for the high accuracy of the "fsolve" # result turns out to be that "fsolve" internally uses a higher value of Digits, and then # rounds its result to the original value of Digits. # # I'll now return Digits to its default value. > Digits:= 10: # # 2. (a) > f:= x -> x^3 + exp(1/x); 3 f := x -> x + exp(1/x) > D(f)(x), (D@@2)(f)(x); 2 exp(1/x) exp(1/x) exp(1/x) 3 x - --------, 6 x + 2 -------- + -------- 2 3 4 x x x # (b) Note that f''(x) > 0 when x > 0. This means that f' is an increasing function. # So for x > 0 there is at most one interval (say 0 < x < a) where f'(x) < 0 (and f is # decreasing) and one (a < x < infinity) where f'(x) > 0 and f is increasing. Therefore # f(x) = c can have at most one solution in the first interval and at most one solution in # the second interval. # # (c) The smallest c is the minimum value of f(x) for x > 0. This will occur at a critical # point of f, i.e. a solution of f'(x) = 0. The fact that f''(x) > 0 for all x > 0 means that if # such a critical point exists, the minimum will occur there. > solve(D(f)(x)=0, x); # That was worth a try, although I didn't really expect it to work. > a:= fsolve(D(f)(x)=0, x); a := .9805089961 > c:= f(a); c := 3.715516997 # (d) Let's first see what "fsolve" by itself gives us. > soln1:= fsolve(f(x)=4, x); soln1 := 1.189115994 # Since this is greater than "a", the other solution should be in the interval 0 .. a. > soln2:= fsolve(f(x)=4, x, x = 0 .. a); soln2 := .8006647233 # (e) To check these solutions, I'll evaluate f at each of them, using a higher setting of # Digits to reduce roundoff error, and compare to 4. > Digits:= 20: > f(soln1) - 4; -9 -.2116296597*10 > f(soln2) - 4; -11 .33701849*10 # Both of the solutions are about as accurate as you could expect with Digits = 10. # In fact, each is a "true" solution rounded to 10 significant digits. > fsolve(f(x)=4, x); 1.1891159940813257727 > fsolve(f(x)=4, x, x = 0 .. a); .80066472330095859085 > Digits:= 10: # # 3 (a) A fixed point is a solution of g(x) = x. > g:= x -> 1.5 * cos(x); g := x -> 1.5 cos(x) > solve(g(x) = x, x); > fp:= fsolve(g(x) = x, x); fp := .9148564784 # So there is a fixed point fp at (approximately) .9148564784. To check whether it's an # attractor or repeller, we check whether |g'(fp)| < 1 or > 1. > D(g)(fp); -1.188712591 # Thus it is a repeller. # # (b) Several starting points could be used. I'll just try one, starting at x = 1. To save space, # I won't show the first 20 iterates, just the 21st to 24th. > X:= 1: for jj from 1 to 20 do X:= evalf(g(X)) od: > for jj from 21 to 24 do X:= evalf(g(X)) od; X := .1252851681 X := 1.488243110 X := .1236892224 X := 1.488540354 # It appears to be settling down to a 2-cycle. # # (c) We can't really prove that most points will approach the 2-cycle, but we can show # that the 2-cycle exists and is an attractor. The points on the 2-cycle will # be solutions of g(g(x)) = x. > x1:= fsolve(g(g(x))=x, x); x1 := 1.488653649 # This was lucky: "fsolve" might have given us the fixed point instead. If so, I could have # asked "fsolve" for a solution in an interval that contained a point of the 2-cycle but not # the fixed point. Now find the other point of the 2-cycle. > x2:= g(x1); x2 := .1230755007 # Now, is the 2-cycle an attractor or repeller? The test is whether |g'(x1) g'(x2)| < 1 or > 1. > D(g)(x1) * D(g)(x2); .2752899284 # So it is an attractor. This means that if we start close to x1 or x2, the iterations will # approach the 2-cycle. # # (d) I'll change the definition of g, and again try an experiment starting at x = 1. This # time no pattern is apparent, even after 100 iterations. I'll just print the 90th to 100th # iterations. Since I want to draw a cobweb diagram, I'll save the results in a table. > g:= x -> 2*cos(x); g := x -> 2 cos(x) > X[0]:= 1: for jj from 1 to 100 do X[jj]:=evalf(g(X[jj-1])) od: > seq(X[jj], jj = 90 .. 100); -.1735617382, 1.969951867, -.7772808092, 1.425646517, .2892813316, 1.916898266, -.6784669592, 1.557071540, .02744871182, 1.999246616, -.8309233368 # The cobweb diagram is made using the "stair" command that was defined in Lesson 4. > stair:= x -> plot([[x,x], [x,g(x)], [g(x), g(x)]]); stair := x -> plot([[x, x], [x, g(x)], [g(x), g(x)]]) > with(plots): # Note that -2 <= g(x) <= 2 for every x, so x = -2 .. 2 is a reasonable interval to use # for the diagram. > p:= plot({x, g(x)}, x = - 2 .. 2): # The first few steps spiral away from the initial point, which was near the fixed point. # But up to the 10th step, there's no evidence that the iterations are approaching anything. > display({p, seq(stair(X[jj]), jj = 0 .. 10) }); # Here are the 80th to 100th steps, still apparently without approaching an attractor. It # seems that this is an instance of chaos. > display({p, seq(stair(X[jj]), jj = 80 .. 100) }); # (e) A way to ensure that a cycle is an attractor is to have g'(x) = 0 for some point on the # cycle. In this case we have g'(0) = 0 (other possibilities are multiples of Pi). To have 0 # on a cycle of period 3, we need (g@@3)(0) = 0. This is an equation which we must # solve for c. I gave "c" a value back in problem 2, so I must "unassign" it to use it as # a symbolic variable. > c:= 'c'; > g:= x -> c * cos(x); c := c g := x -> c cos(x) > eqn:= (g@@3)(0) = 0; eqn := c cos(c cos(c)) = 0 > c1:= fsolve(eqn, c); c1 := 0 # Wrong solution: this makes x=0 a fixed point, not a point on a 3-cycle. We need to give # "fsolve" an interval in which to look. Let's make a graph to see # what interval would be reasonable. > plot(lhs(eqn), c = 0 .. 5); # It looks like the first positive c that works is between 2 and 3. # > c2:= fsolve(eqn, c, c = 2 .. 3); c2 := 2.316112131 >