Chapter 1 The solution of nonlinear Equations∫(x)=0 Consider the physical problem that involves a spherical ball of radius r that is sub- merged to a depth d in water(see Figure 1.1). Assume that the ball is constructed from a variety of longleaf pine that has a density of p=0.638 and that its radius measures r= 10 cm. How much of the ball will be submerged when it is placed in water? The mass Mw of water displaced when a sphere is submerged to a depth d is dx d2(3r-d) and the mass of the ball is Mb= 4rrp/3. Applying Archimedes' law Mo Mb produces the following equation that must be solved 丌(d3-3d2r+4r3p) Figure 1.1 The portion of a sphere of radius r that is to be submerged to depth d
Chapter 1 The Solution of Nonlinear Equations f(x) = 0 Consider the physical problem that involves a spherical ball of radius r that is submerged to a depth d in water(see Figure 1.1). Assume that the ball is constructed from a variety of longleaf pine that has a density of ρ = 0.638 and that its radius measures r = 10 cm. How much of the ball will be submerged when it is placed in water? The mass Mw of water displaced when a sphere is submerged to a depth d is Mw = Z d 0 π(r 2 − (x − r) 2 )dx = πd2 (3r − d) 3 . and the mass of the ball is Mb = 4πr3ρ/3. Applying Archimedes’ law Mw = Mb, produces the following equation that must be solved: π(d 3 − 3d 2 r + 4r 3ρ) 3 = 0. Figure 1.1 The portion of a sphere of radius r that is to be submerged to depth d. 1
1000 Figure 1.2 The cubic y=2552-30d2+d In our case(with r= 10 and p=0.638)this equation becomes x(252-302+=0 The graph of the cubic polynomial y=2552-30d2+ds is shown in Figure 1.2 and from it one can see that the solution lies near the value d=12 The goal of this chapter is to develop a variety of methods for finding numerical approximations for the roots of an equation. For example, the bisection method could be applied to obtain the three roots di=8. 17607212, d2= 11.86150151, and d3 26.31457061. The first root di is not a feasible solution for this problem, because d cannot be negative. The third root ds is larger than the diameter of the sphere and it is not the desired solution. The root d2=11.86150151 lies in the interval [0, 20 and is the proper solution. Its magnitude is reasonable because a little more than one-half of the sphere must be submerge 1.1 Iteration for Solving =g(a A fundamental principle in computer science is iteration. As the name suggests, process is repeated until an answer is achieved. Iterative techniques are used to find roots equations, solutions of linear and nonlinear systems of equations, and solutions of differential equations. In this section we study the process of iteration using repeated substitution A rule or function g(a)for computing successive terms is needed, together with starting value po. Then a sequence of values pk) is obtained using the iterative rule
0 2 4 6 8 10 12 14 16 18 20 −1500 −1000 −500 0 500 1000 1500 2000 2500 3000 d y y=2552−30d2+d3 Figure 1.2 The cubic y = 2552 − 30d 2 + d 3 . In our case (with r = 10 and ρ = 0.638) this equation becomes π(2552 − 30d 2 + d 3 ) 3 = 0. The graph of the cubic polynomial y = 2552 − 30d 2 + d 3 is shown in Figure 1.2 and from it one can see that the solution lies near the value d = 12. The goal of this chapter is to develop a variety of methods for finding numerical approximations for the roots of an equation. For example, the bisection method could be applied to obtain the three roots d1 = −8.17607212, d2 = 11.86150151, and d3 = 26.31457061. The first root d1 is not a feasible solution for this problem, because d cannot be negative. The third root d3 is larger than the diameter of the sphere and it is not the desired solution. The root d2 = 11.86150151 lies in the interval [0, 20] and is the proper solution. Its magnitude is reasonable because a little more than one-half of the sphere must be submerged. 1.1 Iteration for Solving x = g(x) A fundamental principle in computer science is iteration. As the name suggests, a process is repeated until an answer is achieved. Iterative techniques are used to find roots equations, solutions of linear and nonlinear systems of equations, and solutions of differential equations. In this section we study the process of iteration using repeated substitution. A rule or function g(x) for computing successive terms is needed, together with a starting value p0. Then a sequence of values {pk} is obtained using the iterative rule 2
Pk+1= g(pk). The sequence has the pattern starting value) g(po) p2=9(p1) pk=9(pk-1) Pk+1=9(Dk) What can we learn from an unending sequence of numbers? If the numbers tend to a limit, we feel that something has been achieved. But what if the numbers diverge or are periodic? The next example addresses this situation Example 1.1. The iterative rule po= 1 and pk+1= 1.001pk for k=0, 1,..pro- duces a divergent sequence. The first 100 terms look as follows 1.0010=(1.001)(1.000000=1 p2=1.001p1=(1.001)(1.001000)=1.002001, P3=1.001p2=(1.001)(1.002001)=1003003 p100=1.001p-9=(1.001)(1.104012)=1.105116 The process can be continued indefinitely, and it is easily shown that limn-oo Pn =+oo In Chapter 9 we will see that the sequence Ipk is a numerical solution to the differ ential equation y=0.001y. The solution is known to be y(a)=e0.00lzr. Indeed, if we compare the 100th term in the sequence with y(100), we see that P100=1.105116 R 1.105171=c1=y(100) In this section we are concerned with the types of functions g() that produce convergent sequences Pk) 1.1.1 Finding Fixed Points Definition 1.1(Fixed Point). A fired point of a function g(a) is a real number P such that P=g(P) Geometrically, the fixed points of a function y= g(r) are the points of intersectio of y= g(r) and y Definition 1.2(Fixed-point Iteration). The iteration Pn+1=g(pn)for n=0, 1 s called fisced-point iteration
pk+1 = g(pk). The sequence has the pattern p0 (starting value) p1 = g(p0) p2 = g(p1) . . . pk = g(pk−1) pk+1 = g(pk) . . . (1.1) What can we learn from an unending sequence of numbers? If the numbers tend to a limit, we feel that something has been achieved. But what if the numbers diverge or are periodic? The next example addresses this situation. Example 1.1. The iterative rule p0 = 1 and pk+1 = 1.001pk for k = 0, 1, . . . produces a divergent sequence. The first 100 terms look as follows: p1 = 1.001p0 = (1.001)(1.000000) = 1.001000, p2 = 1.001p1 = (1.001)(1.001000) = 1.002001, p3 = 1.001p2 = (1.001)(1.002001) = 1.003003, . . . . . . . . . p100 = 1.001p99 = (1.001)(1.104012) = 1.105116. The process can be continued indefinitely, and it is easily shown that limn→∞ pn = +∞. In Chapter 9 we will see that the sequence {pk} is a numerical solution to the differential equation y 0 = 0.001y. The solution is known to be y(x) = e 0.001x . Indeed, if we compare the 100th term in the sequence with y(100), we see that p100 = 1.105116 ≈ 1.105171 = e 0.1 = y(100). In this section we are concerned with the types of functions g(x) that produce convergent sequences {pk}. 1.1.1 Finding Fixed Points Definition 1.1 (Fixed Point). A fixed point of a function g(x) is a real number P such that P = g(P). Geometrically, the fixed points of a function y = g(x) are the points of intersection of y = g(x) and y = x. Definition 1.2 (Fixed-point Iteration). The iteration pn+1 = g(pn) for n = 0, 1, . . . is called fixed-point iteration. 3
Theorem 1.1. Assume that g is a continuous function and that ipn no is a se- quence generated by fixed-point iteration. If limn-oo pn= P, then P is a fixed point 9(x) Proof. If limn-oo Pn= P, then limn-oo Pn+1= P. It follows from this result, the continuity of 9, and the relation pn+1= g(pn)that g(P)=g(lim pn)= lim g(pn)= lim Pn+1= P (1.2) Therefore, P is a fixed point of g(a) Example 1.2. Consider the convergent iteration po=0.5 and Pk+1 =ePk for k=0, 1, The first 10 terms are obtained by the calculations P1=c-0.50000.606531 0.60031=0.545239 P3=e-0545239=0.579703 0.566409 =0.567560 p10=e-06750=0.566907 The sequence is converging, and further calculations reveal that lim pn2=0.567143. Thus we have found an approximation for the fixed point of the function y =e-r n The following two theorems establish conditions for the existence of a fixed point d the convergence of the fixed- point iteration process to a fixed point Theorem 1.2 Assume that g E Cla, b If the range of the mapping y= g(a)satisfies y E a, b for all x E [a, b],then g has a fixed point in a, b (1.3) Furthermore, suppose that g(a) is defined over (a, b) and that a positive constant K l exists with Ig(a)l< K< l for all E(a, b),then g has a(1.4) unique fixed point P in [a, b
Theorem 1.1. Assume that g is a continuous function and that {pn} ∞ n=0 is a sequence generated by fixed-point iteration. If limn→∞ pn = P, then P is a fixed point of g(x). Proof. If limn→∞ pn = P, then limn→∞ pn+1 = P. It follows from this result, the continuity of g, and the relation pn+1 = g(pn) that g(P) = g( limn→∞ pn) = limn→∞ g(pn) = limn→∞ pn+1 = P. (1.2) Therefore, P is a fixed point of g(x). Example 1.2. Consider the convergent iteration p0 = 0.5 and pk+1 = e −pk for k = 0, 1, . . . . The first 10 terms are obtained by the calculations p1 = e −0.500000 = 0.606531 p2 = e −0.606531 = 0.545239 p3 = e −0.545239 = 0.579703 . . . . . . p9 = e −0.566409 = 0.567560 p10 = e −0.567560 = 0.566907 The sequence is converging, and further calculations reveal that limn→∞ pn = 0.567143 . . . . Thus we have found an approximation for the fixed point of the function y = e −x . The following two theorems establish conditions for the existence of a fixed point and the convergence of the fixed-point iteration process to a fixed point. Theorem 1.2 Assume that g ∈ C[a, b]. If the range of the mapping y = g(x) satisfies y ∈ [a, b] for all x ∈ [a, b], then g has a fixed point in [a, b]. (1.3) Furthermore, suppose that g 0 (x) is defined over (a, b) and that a positive constant K < 1 exists with |g 0 (x)| ≤ K < 1 for all x ∈ (a, b) ,then g has a unique fixed point P in [a, b]. (1.4) 4
Proof of (1.3). If g(a)=a or g(b)=b, the assertion is true. Otherwise, the values of g(a) and g(b)must satisfy g(a E(a, b and g(b)e a, b). The function f(a)= -g(a) has the property that f(a=a-g(a) 1 for all r E a, b, then the iteration Pn=g(pn-1) will not con- verge to P. In this case, P is said to be a repelling fixed point and the iter-(1.7) ation exhibits local divergence
Proof of (1.3). If g(a) = a or g(b) = b, the assertion is true. Otherwise, the values of g(a) and g(b) must satisfy g(a) ∈ (a, b] and g(b) ∈ [a, b). The function f(x) = x − g(x) has the property that f(a) = a − g(a) 0. Now apply Theorem 0.2, the Intermediate Value Theorem, to f(x), with the constant L = 0, and conclude that there exists a number P with P ∈ (a, b) so that f(P) = 0. Therefore, P = g(P) and P is the desired fixed point of g(x). Proof of (1.4). Now we must show that this solution is unique. By way of contradiction, let us make the additional assumption that there exist two fixed points P1 and P2. Now apply Theorem 0.6, the Mean Value Theorem, and conclude that there exists a number d ∈ (a, b) so that g 0 (d) = g(P2) − g(P1) P2 − P1 . (1.5) Next, use the facts that g(P1) = P1 and g(P2) = P2 to simplify the right side of equation (1.5) and obtain g 0 (d) = P2 − P1 P2 − P1 = 1. But this contradicts the hypothesis in (1.4) that|g 0 (x)| 1 for all x ∈ [a, b] ,then the iteration pn = g(pn−1) will not converge to P. In this case, P is said to be a repelling fixed point and the iteration exhibits local divergence. (1.7) 5
Figure 1. 3 The relationship among P, Po, P1, P-Pol and IP-Pl Remark 1. It is assumed that po+P in statement(1.7 Remark 2. Because g is continuous on an interval containing P, it is permissible to use the simpler criterion lg(P)ls Kl in(1.6)and(1, 7), respectively Proof. We first show that the points pn inso all lie in(a, b). Starting with po, we apply Theorem 0.6, the Mean Value Theorem. There exists a value co E(a, b)so that |P-p1|=|9(P)-9(po)=|9(co)(P-p0) g(co川‖P-pl|≤K|P-pol<|P-po Therefore, P1 is no further from P than Po was, and it follows that pi E(a, b)(see Figure 1. 3). In general, suppose that Pn-1 E(a, b);then P-pnl|=|9(P)-9(p1-1)=|g'(cn-1)(P-p2-1) g(cn-1川‖P-pn-1)≤KP-pn-1<|P-p-1) Therefore, Pn E(a, b)and hence, by induction, all the points ipn ingo lie in(a, b) o complete the proof of(1.6), we will show that lim P-pn=0 First, a proof by induction will establish the inequality P-pn|≤K"|P-pol 1.11 The case n= 1 follows from the details in relation(1.8). Using the induction hypothesis IP-Pn-ils kn-IP-pol and the ideas in(1.9), we obtain Pl|≤K|P-p-1≤KK-P-p0|=K|P-pol Thus, by induction, inequality(1.11)holds for all n, Since 0< K< 1, the term K goes to zero as n goes to infinity. Hence ≤l|P-pn|≤linK"|P The limit of IP-pnl is squeezed between zero on the left and zero on the right, so we can conclude that limn-ooP-pnI=0. Thus limn-oo Pn= P and, by Theorem 1.1 the iteration Pn =g(pn-1) converges to the fixed point P. Therefore, statement(1.6) of Theorem 1.3 is proved. We leave statement(1.7) for the reader to investigate
Figure 1.3 The relationship among P, p0, p1, |P − p0| and |P − p1 Remark 1. It is assumed that p0 6= P in statement (1.7) Remark 2. Because g is continuous on an interval containing P, it is permissible to use the simpler criterion |g 0 (P)| ≤ K 1 in (1.6) and (1,7), respectively. Proof. We first show that the points {pn} ∞ n=0 all lie in (a, b). Starting with p0, we apply Theorem 0.6, the Mean Value Theorem. There exists a value c0 ∈ (a, b) so that |P − p1| = |g(P) − g(p0)| = |g 0 (c0)(P − p0)| = |g 0 (c0)||P − p0| ≤ K|P − p0| < |P − p0|. (1.8) Therefore, p1 is no further from P than P0 was, and it follows that p1 ∈ (a, b) (see Figure 1.3). In general, suppose that pn−1 ∈ (a, b); then |P − pn| = |g(P) − g(pn−1)| = |g 0 (cn−1)(P − pn−1)| = |g 0 (cn−1)||P − pn−1)| ≤ K|P − pn−1| < |P − pn−1)|. (1.9) Therefore, pn ∈ (a, b) and hence, by induction, all the points {pn} ∞ n=0 lie in (a, b). To complete the proof of (1.6), we will show that limn→∞ |P − pn| = 0. (1.10) First, a proof by induction will establish the inequality |P − pn| ≤ Kn |P − p0|. (1.11) The case n = 1 follows from the details in relation (1.8). Using the induction hypothesis |P − pn−1| ≤ Kn−1 |P − p0| and the ideas in (1.9), we obtain |P − pn| ≤ K|P − pn−1| ≤ KKn−1 |P − p0| = Kn |P − p0|. Thus, by induction, inequality (1.11) holds for all n, Since 0 < K < 1, the term Kn goes to zero as n goes to infinity. Hence 0 ≤ limn→∞ |P − pn| ≤ limn→∞ Kn |P − p0| = 0. (1.12) The limit of |P − pn| is squeezed between zero on the left and zero on the right, so we can conclude that limn→∞ |P − pn| = 0. Thus limn→∞ pn = P and, by Theorem 1.1, the iteration pn = g(pn−1) converges to the fixed point P. Therefore, statement (1.6) of Theorem 1.3 is proved. We leave statement (1.7) for the reader to investigate. 6
Figure 1.4(a)Monotone convergence when 0< g(P)<1 Figure 1. 4(a) Monotone convergence when -1< g(P)<o Corollary 1.1. Assume that g satisfies the hypothesis given in(1.6)of Theorem 1.3. Bounds for the error involved when using pn to approximate P are given by P-p|≤K"|P- Pol for all n≥1, P-pksA" P-Pol for all n2≥1
Figure 1.4 (a) Monotone convergence when 0 < g0 (P) < 1. Figure 1.4 (a) Monotone convergence when −1 < g0 (P) < 0. Corollary 1.1. Assume that g satisfies the hypothesis given in (1.6) of Theorem 1.3. Bounds for the error involved when using pn to approximate P are given by |P − pn| ≤ Kn |P − p0| for all n ≥ 1, (1.13) and |P − pn| ≤ Kn |P − p0| 1 − K for all n ≥ 1, (1.14) 7
Figure 1.5(a) Monntone divergence when 1<g(P) Figure 1.5 ( b)Divergent oscillation when g(P)<-1 1.1. 2 Graphical Interpretation of Fixed-point Iteration Since we seek a fixed point P to g(ar), it is necessary that the graph of the curve y= g(a) and the line y= a intersect at the point (P, P). Two simple types of convergent iteration, monotone and oscillating, are illustrated in Figure 1.4(a)and(b) respectively To visualize the process, start at po on the x-axis and move vertically to the point (po, P1) g(po)) on the g(a) the point(P1, Pi) on the line y =a. Finally, move vertically downward to Pi on the -axis. The recursion Pn+1 g(pn)is used to construct the point (pn, Pn+1) on the graph, then a horizontal motion locates(Pn+1, Pn+1) on the line y =a, and then a vertical movement ends up at pn+1 on the x-axis. The situation is shown in Figure1. 4
Figure 1.5 (a) Monntone divergence when 1 < g0 (P). Figure 1.5 (b) Divergent oscillation when g 0 (P) < −1. 1.1.2 Graphical Interpretation of Fixed-point Iteration Since we seek a fixed point P to g(x), it is necessary that the graph of the curve y = g(x) and the line y = x intersect at the point (P, P). Two simple types of convergent iteration, monotone and oscillating, are illustrated in Figure 1.4(a) and(b), respectively. To visualize the process, start at p0 on the x -axis and move vertically to the point (p0, p1) = (p0, g(p0)) on the curve y = g(x). Then move horizontally from (p0, p1) to the point (p1, p1) on the line y = x. Finally, move vertically downward to p1 on the x -axis. The recursion pn+1 = g(pn) is used to construct the point (pn, pn+1) on the graph, then a horizontal motion locates (pn+1, pn+1) on the line y = x, and then a vertical movement ends up at pn+1 on the x -axis. The situation is shown in Figure1.4. 8
If lg(P)> 1, then the iteration Pn+1=g(pn)produces a sequence that diverges away from P. The two simple types of divergent iteration, monotone and oscillatin are illustrated in Figure 1.5(a)and(b), respectively Example 1.4. Consider the iteration pn+1= g(pn)when the function g(r)=1+ C-22/ 4 is used. The fixed points can be found by solving the equation x = g(a). The two solutions(fixed points of g) are z=-2 and a=2. The derivative of the function is g(ar)=l-r/2, and there are only two cases to consider Case(i) Case(ii) P=2 Start with P0 =-205 Start with po= 1.6 then get p1=-2.100625 then ge P1=1.96 p2=-220378135 P=1.9996 3=-2.41794414 p3=1.99999996 lim pn lim pn =2. Since lg()l>2 on [-3,1],by The- Since 1g(a)lP and it diverges if we choose po P Example 1.5. Consider the iteration Pn+1 =g(pn) when the function g(a)=2(r- 1)1/2 for z> 1 is used. Only one fixed point P= 2 exists. The derivative is g(a)=l/(r-1)1/2 and g' (2)=1, so Theorem 3. 3 does not apply. There are two cases to consider when the starting value lies to the left or right of P=2 Case(i) Start with po= 1.5 Case(ii Start with po=2.5 then get P1 then get p1=2.44948974 p2=1.28718851 P2=2.40789513 p3=1.07179943 D3=237309514 p4=0.53590832 P4=234358284 p5=2(-0.46409168)12 lim pn=2 nce outside the dor This sequence is converging too slowly g(r), the term ps cannot be comput to the value P= 2; indeed, P1ooo 2.00398714
If |g 0 (P)| > 1, then the iteration pn+1 = g(pn) produces a sequence that diverges away from P. The two simple types of divergent iteration, monotone and oscillating, are illustrated in Figure 1.5 (a) and (b), respectively. Example 1.4. Consider the iteration pn+1 = g(pn) when the function g(x) = 1 + x − x 2/4 is used. The fixed points can be found by solving the equation x = g(x). The two solutions (fixed points of g) are x = −2 and x = 2. The derivative of the function is g 0 (x) = 1 − x/2, and there are only two cases to consider. Case(i) P = −2 Case(ii): P = 2 Start with p0 = −2.05 Start with p0 = 1.6 then get p1 = −2.100625 then get p1 = 1.96 p2 = −2.20378135 p2 = 1.9996 p3 = −2.41794441 p3 = 1.99999996 . . . . . . limn→∞ pn = −∞ limn→∞ pn = 2. Since |g 0 (x)| > 3 2 on [−3, −1], by The- Since |g 0 (x)| P and it diverges if we choose p0 < P. Example 1.5. Consider the iteration pn+1 = g(pn) when the function g(x) = 2(x − 1)1/2 for x ≥ 1 is used. Only one fixed point P = 2 exists. The derivative is g 0 (x) = 1/(x − 1)1/2 and g 0 (2) = 1, so Theorem 3.3 does not apply. There are two cases to consider when the starting value lies to the left or right of P = 2. Case(i) Start with p0 = 1.5 Case(ii): Start with p0 = 2.5 then get p1 = 1.41421356 then get p1 = 2.44948974 p2 = 1.28718851 p2 = 2.40789513 p3 = 1.07179943 p3 = 2.37309514 p4 = 0.53590832 p4 = 2.34358284 . . . . . . p5 = 2(−0.46409168)1/2 . limn→∞ pn = 2. Since p4 lies outside the domain of This sequence is converging too slowly g(x), the term p5 cannot be computed to the value P = 2; indeed, P1000 = 2.00398714. 9
1.1.3 Absolute and relative Error Considerations In Example 1.5, case(ii), the sequence converges slowly, and after 1000 iterations the three consecutive terms are F1000=2.00398714,P0o1=2.00398317,and1002=2.00397921 This should not be disturbing; after all, we could compute a few thousand more terms and find a better approximation! But what about a criterion for stopping the iteration? Notice that if we use the difference between consecutive terms p101-p100=|2.00398317-2.003979211=0.0000396 Yet the absolute error in the approximation P1ooo is known to be |P-po|=2.000002.00398714=0.00398714 This is about 1000 times larger than P1001-P1002 and it shows that closeness of consecutive terms does not guarantee that accuracy has been achieved. But it is usually the only criterion available and is often used to terminate an iterative procedure Program 1.1(Fixed-Point Iteration). To approximate a solution to the equation =g() starting with the initial guess Po and iterating Pn+1= g(pn Function[k, p, err, P]=fixpt(g, po, tol, max1) %%%%%%%% %Input- g is the iteration function input as a string g po is the initial guess for the fixed tol is the tolerance max1 is the maximum number of iterations %Output- k is the number of iterations that were carried out -p is the approximation to the fixed point err is the error in the approximation P contains the sequence ( pn) P(1)=po; for k=2. max1 P(k)=feval(g, P(k-1)) err=abs(P(k)-P(k-1 relerr=err/abs(P(k))+eps) if(err<tol)I(relerr<tol), break;end end if k== max1 disp(maximum number of iterations exceeded) d P=P Remark. When using the user-defined function fixpt, it is necessary to input the M-file g m as a string: 'g'(see MATLAB Appendix) 10
1.1.3 Absolute and Relative Error Considerations In Example 1.5, case (ii), the sequence converges slowly, and after 1000 iterations the three consecutive terms are P1000 = 2.00398714, P1001 = 2.00398317, and P1002 = 2.00397921. This should not be disturbing; after all, we could compute a few thousand more terms and find a better approximation! But what about a criterion for stopping the iteration? Notice that if we use the difference between consecutive terms, |p1001 − p1002| = |2.00398317 − 2.00397921| = 0.00000396. Yet the absolute error in the approximation P1000 is known to be |P − p1000| = |2.00000000 − 2.00398714| = 0.00398714. This is about 1000 times larger than |p1001 − p1002| and it shows that closeness of consecutive terms does not guarantee that accuracy has been achieved. But it is usually the only criterion available and is often used to terminate an iterative procedure. Program 1.1 (Fixed-Point Iteration). To approximate a solution to the equation x = g(x) starting with the initial guess p0 and iterating pn+1 = g(pn). Function [k, p, err, P] =fixpt(g, po, tol, max1) % Input − g is the iteration function input as a string 0 g 0 % − po is the initial guess for the fixed point % − tol is the tolerance % − max1 is the maximum number of iterations %Output− k is the number of iterations that were carried out % − p is the approximation to the fixed point % − err is the error in the approximation % − P contains the sequence {pn} P(1)= po; for k=2:max1 P(k)=feval(g, P(k−1)); err=abs(P(k)−P(k−1)); relerr=err/(abs(P(k))+eps); p=P(k); if (err<tol) | (relerr<tol), break; end end if k == max1 disp(’maximum number of iterations exceeded’) end P=P’; Remark. When using the user-defined function fixpt, it is necessary to input the M-file g.m as a string: ’g’ (see MATLAB Appendix). 10
