Page 286 - MATLAB an introduction with applications

Optimization ——— 271

>> x0=[0 0];TolX=1e–4;TolFun=1e–6;MaxIter=100;
>> [x0,g0, xx]=newtons(gn,x0,TolX,MaxIter)
x0 =
3.6667 2.3333
g0 =
1.0e–015 *
0 –0.4441
xx =
3.6667 2.3333
3.6667 2.3333

function g= jacob(f,x,h,varargin)
%Jacobian of f(x)
if nargin<3, h=.0001; end
N= length(x); h2= 2*h; %h12=12*h;
x=x(:).’; I= eye(N);
for n=1:N
f1=feval(f,x+I(n,:)*h,varargin{:});
f2=feval(f,x–I(n,:)*h,varargin{:});
f3=feval(f,x+I(n,:)*h2,varargin{:});
f4=feval(f,x–I(n,:)*h2,varargin{:});
f12=(f1–f2)/h2;
f12=(8*(f1–f2)–f3+f4)/h12;
g(:,n)=f12(:);
end
if sum(sum(isnan(g)))==0&rank(g)<N
format short e
fprintf(‘At x=%12.6e, Jacobian singular with J=’,x);
disp(g); format short;
end

function [x,fx,xx]=newtons(f,x0,TolX,MaxIter,varargin)

% newtons.m to solve a set of nonlinear eqs
% input:f=a 1st-order vector ftn equivalent to a set of equations
% x0=the initial guess of the solution
% TolX=the upper limit of |x(k)–x(k–1)|
% MaxIter=the maximum # of iteration
% Output: x=the point which the algorithm has reached
% fx=f(x(last))
% xx=the history of x