# HW8 - Homework 8 Solution hw8_solution.m clear all tol =...

This preview shows page 1. Sign up to view the full content.

This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: Homework 8 Solution hw8_solution.m clear all tol = 1e‐20; xGuess = 1.2; yGuess = .7; outputVec = hw8delG_ Solution([xGuess,yGuess]); dfDx = outputVec(1); dfDy = outputVec(2); err = dfDx^2 + dfDy^2; xGuessStore = [xGuess]; yGuessStore = [yGuess]; while(err > tol) % We need to take a step. So we use optimize our % single‐dimensional function to find the value of a. % Don't use [‐2, 2] for bounds, because we don't know what a should be newA = goldenSectionMax_ Solution(@(a) hw8singleParamG_ Solution(a,xGuess,yGuess,dfDx,dfDy),... ‐100,100); %plot([‐100:100],hw8singleParamG_ Solution([‐100:100],xGuess,yGuess,dfDx,dfDy)) xGuess = xGuess + dfDx*newA; yGuess = yGuess + dfDy*newA; outputVec = hw8delG_ Solution([xGuess,yGuess]); dfDx = outputVec(1); dfDy = outputVec(2); err = sqrt(dfDx^2 + dfDy^2); xGuessStore = [xGuessStore xGuess]; yGuessStore = [yGuessStore yGuess]; end solution = [xGuess yGuess] fVal = hw7funcG_ Solution(xGuess,yGuess) firstIterations = [xGuessStore(1:7);yGuessStore(1:7)] % PART 3: nrRes = newtonRaphson_ Solution(@hw8delG_ Solution,@hw8hessianG_ Solution,[1;1]) % PART 4: % We can play around with the function in 2‐d to make sure we can get it % right based on that simplification. Note we would then need to change % the gradient and Hessians as well. I've not included that here. initGuess = [0;0;1]; nrRes = newtonRaphson_ Solution(@hw8delW_ Solution,@hw8HessianW_ Solution,initGuess) % We can also make sure that Matlab agrees with us: % (Note Matlab looks for a minimum, so we use our @(x) formalism to change % our function so that we're looking for a minimum fminsearch(@(x) ‐1*hw8functionW_ Solution(x),initGuess) function outputVec = hw8delG_ Solution(inputVec) x = inputVec(1); y = inputVec(2); dfDx = ‐8*x + 2*y; dfDy = 2*x ‐ 2*y; outputVec = [dfDx; dfDy]; function f = hw8funcG_ Solution(x,y) % x = inputVec(1); % y = inputVec(2); f = ‐4*x.^2 + 2*x.*y ‐ y.^2; function maxVal = goldenSectionMax_ Solution(optFunc,x_l,x_u) % Golden section method % Find the function values at our upper and lower bounds (not actually % necessary) x_uFval =optFunc(x_u); x_lFval = optFunc(x_l); % The Golden ratio phi = (sqrt(5)‐1)/2; % Our starting interval width, used to define d at this step intervalWidth = x_u ‐ x_l; x_1 = x_l + phi*intervalWidth; x_2 = x_u ‐ phi*intervalWidth; % f(x_1) and f(x_2) used to make decisions for golden section x_1fval = optFunc(x_1); x_2fval = optFunc(x_2); myTol = 1e‐6; while (intervalWidth > myTol) if (x_2fval > x_1fval) % Eliminate top portion of interval x_u = x_1; x_uFval = x_1fval; % Make x2 be the new x_1 x_1 = x_2; x_1fval = x_2fval; % Set a new x_2 intervalWidth = x_u ‐ x_l; x_2 = x_u ‐ phi*intervalWidth; x_2fval = optFunc(x_2); else % Eliminate the bottom portion of the interval x_l = x_2; x_lFval = x_2fval; % Make x1 be the new x2 x_2 = x_1; x_2fval = x_1fval; % Set a new x_1 intervalWidth = x_u ‐ x_l; x_1 = x_l + phi*intervalWidth; x_1fval = optFunc(x_1); end end % Our best estimate is probably the middle of this interval maxVal = (x_u + x_l)/2; function ans = hw8singleParamG_ Solution(a,xi,yi,dfiDx,dfiDy) ans = ‐4*(xi+dfiDx*a).^2 + 2*(xi+dfiDx*a).*(yi+dfiDy*a)‐(yi+dfiDy*a).^2; function outputVec = hw8delG_ Solution(inputVec) x = inputVec(1); y = inputVec(2); dfDx = ‐8*x + 2*y; dfDy = 2*x ‐ 2*y; outputVec = [dfDx; dfDy]; Newton Raphson function root = newtonRaphson_ Solution(rootFun,jacFun,initGuess) % Set our solution tolerance and information about the initial guess tol = 1e‐10; guess = initGuess; guessFval = rootFun(guess); % Iterate while the guess is not within tolerance of zero while(max(abs(guessFval) > tol)) % Use the Jacobian and Matlab's built‐in GE methods to solve our linear % system jacVal = jacFun(guess); newGuess = jacVal\(jacVal*guess ‐ guessFval); newGuessFval = rootFun(newGuess); % Replace the old guess with the new guess guess = newGuess; guessFval = newGuessFval; end % What remains is the solution. root = guess; function H = hw8hessianG_ Solution(x,y) H = [‐8 2; 2 ‐2]; function outVec = hw8delW_ Solution(inVec) outVec = zeros(3,1); % outVec = zeros(2,1); x1 = inVec(1); x2 = inVec(2); x3 = inVec(3); outVec(1) = ‐sin(x1+x2) + cos(x1); outVec(2) = ‐sin(x1+x2) ‐ sin(x2); outVec(3) = ‐0.2*x3.*exp(‐01.*x3.^2); function fVal = hw8functionW_ Solution(inputVec) x1 = inputVec(1); x2 = inputVec(2); x3 = inputVec(3); fVal = cos(x1+x2) + sin(x1) + cos(x2) + exp(‐0.1*x3.^2); function H = hw8HessianW_ Solution(inputVec) x1 = inputVec(1); x2 = inputVec(2); x3 = inputVec(3); H = zeros(3); H(1,1) = ‐cos(x1+x2)‐sin(x1); H(1,2) = ‐cos(x1+x2); H(1,3) = 0; H(2,1) = ‐cos(x1+x2); H(2,2) = ‐cos(x1+x2)‐cos(x2); H(2,3) = 0; H(3,1) = 0; H(3,2) = 0; H(3,3) = (‐0.2*x3).^2.*exp(‐0.1*x3.^2) ‐ 0.2*exp(‐0.1*x3.^2); % H = zeros(2); % H(1,1) = ‐cos(x1+x2)‐sin(x1); % H(1,2) = ‐cos(x1+x2); % H(2,1) = ‐cos(x1+x2); % H(2,2) = ‐cos(x1+x2)‐cos(x2); function fVal = hw8functionWforMesh_ Solution(x1,x2,x3) fVal = cos(x1+x2) + sin(x1) + cos(x2) + exp(‐0.1*x3.^2); ...
View Full Document

## This note was uploaded on 05/01/2011 for the course CHBE 2120 taught by Professor Gallivan during the Spring '07 term at Georgia Tech.

Ask a homework question - tutors are online