Ma 449: Numerical Applied Mathematics Model Solutions to Homework Assignment 10 Prof. Wickerhauser Exercise 1* of Section 8.1, pp.427--428. (a) f'(x) = 6x^2-18x+12 = 6(x-1)(x-2), so f(x) increases on (-infty, 1) where f'(x) > 0; decreases on (1, 2) where f'(x) < 0; increases on (2, infty) where f'(x) > 0. (b) f'(x)=1/(x+1)^2 >0, so f(x) increases on R - {-1} (c) This is the reciprocal of the function in part (b). It decreases where the part (b) function increases, so f(x) decreases on R - {-1}. (d) f(x)=x^x = exp(x log x), so f'(x) = (1+log x) exp(x log x), so f(x) decreases on (0, 1/e), where f'(x)<0 increases on (1/e, infty), where f'(x)>0. Exercise 2(b) of Section 8.1, pp.427--428. Observe that cos(x) is decreasing on [0, pi] and increasing on [pi, 2*pi], and that 0 < pi < 4 < 2*pi. Conclude that f is unimodal on [0,4] with a unique minimum at pi. Exercise 7(d*) of Section 8.1, pp.427--428. C PROGRAM: /* Golden ratio search with N iterations for the minimum of a function f(x) unimodal on [a0,b0] */ #include #include #include double f(double x) { return -sin(x)-x+x*x/2; } double a0=0.8, b0=1.6; int N=20; /* number of iterations */ main() { double a, b, c, d, fa, fb, fc, fd; const double r=(sqrt(5.0)-1)/2; /* golden ratio */ int k; a=a0; b=b0; fa=f(a); fb=f(b); /* initialize */ c=r*a+(1-r)*b; d=r*b+(1-r)*a; fc=f(c); fd=f(d); /* a < c < d < b */ for(k=0; k fd ) { /* ==> min lies in [c,b] */ a=c; c=d; d=r*b+(1-r)*a; /* b is still the same */ fa=fc; fc=fd; fd=f(d); /* f(b) is still the same */ } else { /* min lies in [a,d] */ b=d; d=c; c=r*a+(1-r)*b; /* a is still the same */ fb=fd; fd=fc; fc=f(c); /* f(a) is still the same */ } } return 0; } OUTPUT: k=0: a[0]=0.8000, b[0]=1.6000 k=1: a[1]=1.1056, b[1]=1.6000 k=2: a[2]=1.1056, b[2]=1.4111 k=3: a[3]=1.2223, b[3]=1.4111 k=4: a[4]=1.2223, b[4]=1.3390 k=5: a[5]=1.2669, b[5]=1.3390 k=6: a[6]=1.2669, b[6]=1.3115 k=7: a[7]=1.2669, b[7]=1.2944 k=8: a[8]=1.2774, b[8]=1.2944 k=9: a[9]=1.2774, b[9]=1.2879 k=10: a[10]=1.2814, b[10]=1.2879 k=11: a[11]=1.2814, b[11]=1.2854 k=12: a[12]=1.2830, b[12]=1.2854 k=13: a[13]=1.2830, b[13]=1.2845 k=14: a[14]=1.2830, b[14]=1.2839 k=15: a[15]=1.2833, b[15]=1.2839 k=16: a[16]=1.2833, b[16]=1.2837 k=17: a[17]=1.2833, b[17]=1.2835 k=18: a[18]=1.2834, b[18]=1.2835 k=19: a[19]=1.2834, b[19]=1.2835 Exercise 8(d) of Section 8.1, pp.427--428. Modify the Fibonacci method Matlab function from the textbook's web site: fib.m from http://math.fullerton.edu/mathews/software/NumericalMethods7.1.zip function y=fib(n) %Generates fibonacci numbers for the program fibonacci fz(1)=1;fz(2)=1; for k=3:n fz(k)=fz(k-1)+fz(k-2); end y=fz(n); fibonacci10.m (modified from fibonacci.m) from http://math.fullerton.edu/mathews/software/NumericalMethods7.1.zip function X=fibonacci10(f,a,b,tol,e) %Input - f, the object function % - a, the left endpoint of the interval % - b, the right endpoint of the interval % - tol, length of uncertainty % - e, distinguishability constant %Output - X, x and y coordinates of minimum %Note this function calls the m-file fib.m %If f is defined as an M-file function use the @ notation % call X=fibonacci(@f,a,b,tol,e). %If f is defined as an anonymous function use the % call X=fibonacci(f,a,b,tol,e). % NUMERICAL METHODS: Matlab Programs % (c) 2004 by John H. Mathews and Kurtis D. Fink % Complementary Software to accompany the textbook: % NUMERICAL METHODS: Using Matlab, Fourth Edition % ISBN: 0-13-065248-2 % Prentice-Hall Pub. Inc. % One Lake Street % Upper Saddle River, NJ 07458 %%%%%%%%% CHANGE OLD: %%Determine n %i=1; %F=1; %while F<=(b-a)/tol % F=fib(i); % i=i+1; %end %%%%%%%%%%% ...TO NEW: % Assume n=10 as the problem requests: i=10+1; %%%%%%%%%%%%% %Initialize values n=i-1; A=zeros(1,n-2);B=zeros(1,n-2); A(1)=a; B(1)=b; c=A(1)+(fib(n-2)/fib(n))*(B(1)-A(1)); d=A(1)+(fib(n-1)/fib(n))*(B(1)-A(1)); k=1; %Compute Iterates while k<=n-3 if f(c)>f(d) A(k+1)=c; B(k+1)=B(k); c=d; d=A(k+1)+(fib(n-k-1)/fib(n-k))*(B(k+1)-A(k+1)); else A(k+1)=A(k); B(k+1)=d; d=c; c=A(k+1)+(fib(n-k-2)/fib(n-k))*(B(k+1)-A(k+1)); end k=k+1; end %Last iteration using distinguishability constant e if f(c)>f(d) A(n-2)=c; B(n-2)=B(n-3); c=d; d=A(n-2)+(0.5+e)*(B(n-2)-A(n-2)); else A(n-2)=A(n-3); B(n-2)=d; d=c; c=A(n-2)+(0.5-e)*(B(n-2)-A(n-2)); end %Output: Use midpoint of last interval for abscissa if f(c)>f(d) a=c;b=B(n-2); else a=A(n-2);b=d; end X=[(a+b)/2 f((a+b)/2)]; Implement the function of section 8.1, exercise 7(d), p.428: function y=f(x) y = 23 + x.*x.*(x-5); OUTPUT: >> fibonacci10(@f,1,5,.00001,.000001) ans = 3.3636 4.4861 Algorithm 2(d) of Section 8.1, p.430. This is a straightforward application of the Fibonnaci search of Progeam 8.2, p.423, on the same function f() as in Exercise 8 of Section 8.1, pp.427--428 above. OUTPUT: >> fibonacci(@f,1,5,.00001,.000001) ans = 3.3333 4.4815 Algorithm 3(d) of Section 8.1, p.430. Copy Program 8.3 on p.425, from the textbook web site: quadmin.m, from http://math.fullerton.edu/mathews/software/NumericalMethods7.1.zip function [p,yp,dp,dy,P] = quadmin(f,a,b,delta,epsilon) %Input - f is the object function % - a and b are the endpoints of the interval % - delta is the tolerance for the abscissas % - epsilon is the tolerance for the ordinates %Output - p is the abscissa of the minimum % - yp is the ordinate of the minimum % - dp is the error bound for p % - dy is the error bound for yp % - P is the vector of iterations %If f is defined as an M-file function use the @ notation % call [p,yp,dp,dy,P] = quadmin(@f,a,b,delta,epsilon). %If f is defined as an anonymous function use the % call [p,yp,dp,dy,P] = quadmin(f,a,b,delta,epsilon). % NUMERICAL METHODS: Matlab Programs % (c) 2004 by John H. Mathews and Kurtis D. Fink % Complementary Software to accompany the textbook: % NUMERICAL METHODS: Using Matlab, Fourth Edition % ISBN: 0-13-065248-2 % Prentice-Hall Pub. Inc. % One Lake Street % Upper Saddle River, NJ 07458 p0 = a; maxj = 20; maxk = 30; big = 1e6; err = 1; k = 1; P(k) = p0; cond = 0; h = 1; if (abs(p0)>1e4), h = abs(p0)/1e4; end while (kepsilon & cond~=5) f1 = (f(p0+0.00001)-f(p0-0.00001))/0.00002; if (f1>0), h = -abs(h); end p1 = p0 + h; p2 = p0 + 2*h; pmin = p0; y0 = f(p0); y1 = f(p1); y2 =f(p2); ymin = y0; cond = 0; j = 0; %Determine h so that y1delta & cond==0) if (y0<=y1), p2 = p1; y2 = y1; h = h/2; p1 = p0 + h; y1 = f(p1); else if (y2big | abs(p0)>big), cond=5; end end if (cond==5), pmin = p1; ymin =f(p1); else %Quadratic interpolation to find yp d = 4*y1-2*y0-2*y2; if (d<0), hmin = h*(4*y1-3*y0-y2)/d; else hmin = h/3; cond = 4; end pmin = p0 + hmin; ymin =f(pmin); h = abs(h); h0 = abs(hmin); h1 = abs(hmin-h); h2 = abs(hmin-2*h); %Determine magnitude of next h if (h0big | abs(pmin)>big), cond=5; end %Termination test for minimization e0 = abs(y0-ymin); e1 = abs(y1-ymin); e2 = abs(y2-ymin); if (e0~=0 & e0> quadmin(@f,1,5,.00001,.00001) ans = 3.3333 Exercise 2 of Section 8.2, pp.444-445. MatLab commands: >> b=[2 3]; g=[1 1]; w=[5 2]; >> m=(b+g)/2; r=2*m-w; e=2*r-m; >> xy=[b; g; w; b; r; g; e; b; m; r; e]; >> x=xy(:,1); y=xy(:,2); >> plot(x,y) See the plot at ./8-2-2.jpg Exercise 6 of Section 8.2, pp.444-445. % Powell's method -- minimization u1=[1 0]; u2=[0 1]; % Initial unit coordinate vectors p0=[1/2 1/3]; % Initial point f(x,y)=x^3+y^3-3x-3y+5 % Compute p1 g1(t)=f(p0+t*u1)= (1/2+t)^3+(1/3)^3-3(1/2+t)-3(1/3)+5 g1'(t) = 3(1/2+t)^2-3 = 0 ==> t=1/2 (argmin) or t= -3/2 (argmax) t1 = 1/2 % since we want the minimum, as in example 8.7. p1 = p0+t1*u1 = [1 1/3] % Compute p2 g2(t)=f(p1+t*u2)= (1)^3+(1/3+t)^3-3(1)-3(1/3+t)+5 g2'(t) = 3(1/3+t)^2-3 = 0 ==> t=2/3 (argmin) or t= -4/3 (argmax) t2 = 2/3 % since we want the minimum, as in example 8.7. p2 = p1+t2*u2 = [1 1] x1 = p2 = [1 1] % Answer Algorithm 2(a) of Section 8.2, p.445. Note: Program 8.4 on p.441 of the text is available online from http://math.fullerton.edu/mathews/software/NumericalMethods7.1.zip However, it has two bugs and does not work correctly with the current version of Matlab. Use instead the sample "fixed" nelder.m from the class website: http://www.math.wustl.edu/~victor/classes/ma449/nelder.m %Call [V0,y0,dV,dy]=nelder(@F,V,min1,max1,epsilon,show) % Corrected by MVW, 2007-11-14 % %Input - F is input as an M-file function % - V is a 4x3 matrix containing the starting simplex % CORRECTED % - min1 & max1 are minimum and maximum number of iterations % - epsilon is the tolerance % - show == 1 displays iterations (P and Q) %Output - V0 is the vertex forthe minimum % - y0 is the function value F(V0) % - dV is the size of the final simplex % - dy is the error bound for the minimum % - P is a matrix containing the vertex iterations % - Q is an array containing the iterations for F(P) ... Create the file "f.m" with the following function definition: function out = f(in) x = in(1); y=in(2); z=in(3); out = 2*x.^2 +2*y^2 + z^2 - 2*x*y + y*z - 7*y - 4*z; Start Matlab and run the following commands: V = [1 1 1; 0 1 0; 1 0 1; 0 0 1]; % note: transpose to suit fixed nelder() mini = 10; maxi = 100000; % limits for number of iterations epsilon = 0.5e-8; % tolerance, or diameter of target simplex show = 0; % do not display iterations format('long') % display 15 decimal places of output [V0,y0,dV,dy]=nelder(@f,V,mini,maxi,epsilon,show) OUTPUT: V0 = 1.000010421530515 1.999968981847138 0.999998570523933 y0 = -8.999999997165638 dV = 5.403293974892698e-05 dy = 3.305942186671018e-09 Exercise 1(b*) of Section 8.3, p.456 Hint: Use Macsyma. The commands and their outputs are: f(x,y):= 100*(y-x^2)^2 + (1-x)^2; f(1/2,4/3); 2117 ---- 18 diff(f(x,y),x,1); 2 - 400 x (y - x ) - 2 (1 - x) at(%,[x=1/2,y=4/3]); 653 - --- 3 diff(f(x,y),y,1); 2 200 (y - x ) at(%,[x=1/2,y=4/3]); 650 --- 3 Thus: Grad f(x,y) = ( 200*(-2*x)*(y-x^2)-2*(1-x), 200*(y-x^2) ) Grad f(1/2, 4/3) = ( -653/3, 650/3 ) ~ (-217.6667, 216.6667) Exercise 2(b) of Section 8.3, p.456 Search direction at (1/2,4/3) is Grad f(1/2, 4/3) = ( -653/3, 650/3 ) -- normalized to unit length, this is ~( -0.7087, 0.7055) Find P_1 along the line through (1/2,4/3) in the search direction: L(t) = (1/2,4/3) + t*( -653/3, 650/3 ) = ( 1/2-653*t/3, 4/3-650*t/3 ) f(L(t)) = Using Macsyma, the commands and their outputs are x(t):= 1/2 + t*(-653/3); y(t) := 4/3 + t*(650/3); g(t) := f(x(t),y(t)); diff(g(t),t,1); 653 t 1 1 653 t 1306 (----- + -) 1306 (- - -----) 3 2 2 3 650 ---------------- + 200 (---------------- + ---) 3 3 3 650 t 1 653 t 2 4 (----- - (- - -----) + -) 3 2 3 3 Expanding this as a polynomial gives expand(%); 3 2 72729854112400 t 111122185400 t 466912154 t 848909 ----------------- - --------------- + ----------- + ------ 81 9 27 9 Get floating-point approximations of the coefficients with float(%); 8.978994334864198E+11 * t^3 - 1.234690948888889E+10 * t^2 + 1.7293042740740744E+7 * t + 94323.22222222222 Find the unique real root t0 with Newton's method as implemented in Macsyma: load(newton1); /usr/share/maxima/5.16.3/share/numeric/newton1.mac newton(%,t,0.5,1e-12); .01119116736667045 Finally, evaluate the line at this value t0 of t: x(.01119116736667045); - 1.935944096811934 y(.01119116736667045); 3.758086262778597 g(.01119116736667045); 8.630185446383528 This gives P1 ~ ( - 1.93594, 3.75808 ), where the value of f(x(t0),y(t0)) = g(t0) is about 8.63. Exercise 3(b*) of Section 8.3, p.456 Hint: Use Macsyma. The commands and their outputs are: diff(f(x,y),x,2); 2 2 - 400 (y - x ) + 800 x + 2 at(%,[x=1/2,y=4/3]); 694 - --- 3 diff(f(x,y),x,1, y,1); - 400 x at(%,[x=1/2,y=4/3]); - 200 diff(f(x,y),y,1, x,1); - 400 x at(%,[x=1/2,y=4/3]); - 200 diff(f(x,y),y,2); 200 at(%,[x=1/2,y=4/3]); 200 The Hessian matrix at a general point (x,y) is thus Hf(x,y) = d^2f/dx^2 d^2f/dxdy = -400(y-x^2)+800x^2+2 -400x d^2f/dydx d^2f/dy^2 -400x 200 Hence Hf(1/2,4/3) = -694/3 -200 -200 200 Exercise 4(b) of Section 8.3, p.456 The second-degree Taylor polynomial is expressed in terms of the gradient and Hessian as follows: f(x-1/2,y-4/3) = f(1/2, 4/3) + (x,y)*grad f(1/2,3/4) + (1/2)*(x,y)*Hf(1/2,4/3)*(x,y)' Using the results from Exercises 1(b) and 3(b) in this section, we get f(x-1/2,y-4/3) = (2117/18) + (-653/3)*x + (650/3)*y + (-694/6)*x^2 + (-200)*x*y + (100)*y^2] ~ 117.61 - 217.67 x + 216.66 y - 115.67 x^2 - 200 x y + 100 y^2 Exercise 6(b) of Section 8.3, p.456 The direction vector for the modified Newton's method is s = - Hf(1/2,4/3)^(-1) * grad f(1/2,4/3) Evaluate this using the values computed in Macsyma above gives Hf:matrix([-694/3, -200], [ -200, 200]); [ 694 ] [ - --- - 200 ] [ 3 ] [ ] [ - 200 200 ] invert(Hf); [ 3 3 ] [ - ---- - ---- ] [ 1294 1294 ] [ ] [ 3 347 ] [ - ---- ------ ] [ 1294 129400 ] gradf:matrix([-653/3], [650/3]); [ 653 ] [ - --- ] [ 3 ] [ ] [ 650 ] [ --- ] [ 3 ] Matrix multiplication in Macsyma is denoted with a period: invert(Hf).gradf; [ 3 ] [ ---- ] [ 1294 ] [ ] [ 8429 ] [ ---- ] [ 7764 ] Floating-point approximations to the direction vector coefficients: float(%); [ 0.00231839258114374 ] [ ] [ 1.085651725914477 ] The actual step in Newton's method is the negative of this vector, but in the modified method the sign is not relevant. This modified Newton search vector should be compared with the search vector ( -0.7087,0.7055) from the gradient method. For this particular example, they are quite different.