VL14: Iterative Verfahren zur Nullstellenbestimmung

Contents

Bisektion zur Berechnung von sqrt(2)

clear all;
close all;
f = @(x) x.^2 - 2;


a = 1;
b = 2;
k = 0;
while abs(a-b)>eps
    x = (a+b)/2;
    if sign(f(x)) == sign(f(b))
        b = x;
    else
        a = x;
    end
    k = k+1;
end
fprintf('Nullstelle: %1.8f\n',x)
fprintf('Anzahl der Iterationen: %d\n',k)

Nullstelle: 1.41421356
Anzahl der Iterationen: 52

Newton-Iteration

fprime = @(x) 2.*x;

k=0;
x=1.5;
xprev = 1;
while abs(x - xprev) > eps
    xprev = x;
    x = x - f(x) / fprime(x);
    k = k+1;
end
fprintf('Nullstelle: %1.8f\n',x)
fprintf('Anzahl der Iterationen: %d\n',k)

Nullstelle: 1.41421356
Anzahl der Iterationen: 5

Sekantenverfahren

a = 1;
b = 2;
k=0;
while abs(b-a) > eps
    c = a;
    a = b;
    b = b + (b-c) / (f(c) / f(b) -1);
    k = k+1;
end
fprintf('Nullstelle: %1.8f\n',a)
fprintf('Anzahl der Iterationen: %d\n',k)

Nullstelle: 1.41421356
Anzahl der Iterationen: 7

Inverse quadratische Interpolation

a = 1;
b=1.5;
c=2;

k=0;
while abs(c-b) > eps
    x = polyinterp([f(a),f(b),f(c)],[a,b,c],0);
    a = b;
    b = c;
    c = x;
    k = k+1;
end
fprintf('Nullstelle: %1.8f\n',a)
fprintf('Anzahl der Iterationen: %d\n',k)

Nullstelle: 1.41421356
Anzahl der Iterationen: 6

Matlab-Implementation von fzero:

x = fzerotx(f,[1,2])

Number of Iterations in fzerotx: 7

x =

    1.4142

Published with MATLAB® R2014a

function b = fzerotx(F,ab,varargin)
%FZEROTX  Textbook version of FZERO.
%   x = fzerotx(F,[a,b]) tries to find a zero of F(x) between a and b.
%   F(a) and F(b) must have opposite signs.  fzerotx returns one
%   end point of a small subinterval of [a,b] where F changes sign.
%   Arguments beyond the first two, fzerotx(F,[a,b],p1,p2,...),
%   are passed on, F(x,p1,p2,..).
%
%   Examples:
%      fzerotx(@sin,[1,4])
%      F = @(x) sin(x); fzerotx(F,[1,4])

% Initialize.
a = ab(1);
b = ab(2);
fa = F(a,varargin{:});
fb = F(b,varargin{:});
if sign(fa) == sign(fb)
   error('Function must change sign on the interval')
end
c = a;
fc = fa;
d = b - c;
e = d;

k=0;
% Main loop, exit from middle of the loop
while fb ~= 0

   % The three current points, a, b, and c, satisfy:
   %    f(x) changes sign between a and b.
   %    abs(f(b)) <= abs(f(a)).
   %    c = previous b, so c might = a.
   % The next point is chosen from
   %    Bisection point, (a+b)/2.
   %    Secant point determined by b and c.
   %    Inverse quadratic interpolation point determined
   %    by a, b, and c if they are distinct.

   if sign(fa) == sign(fb)
      a = c;  fa = fc;
      d = b - c;  e = d;
   end
   if abs(fa) < abs(fb)
      c = b;    b = a;    a = c;
      fc = fb;  fb = fa;  fa = fc;
   end

   % Convergence test and possible exit
   m = 0.5*(a - b);
   tol = 2.0*eps*max(abs(b),1.0);
   if (abs(m) <= tol) | (fb == 0.0)
      break
   end

   % Choose bisection or interpolation
   if (abs(e) < tol) | (abs(fc) <= abs(fb))
      % Bisection
      d = m;
      e = m;
   else
      % Interpolation
      s = fb/fc;
      if (a == c)
         % Linear interpolation (secant)
         p = 2.0*m*s;
         q = 1.0 - s;
      else
         % Inverse quadratic interpolation
         q = fc/fa;
         r = fb/fa;
         p = s*(2.0*m*q*(q - r) - (b - c)*(r - 1.0));
         q = (q - 1.0)*(r - 1.0)*(s - 1.0);
      end;
      if p > 0, q = -q; else p = -p; end;
      % Is interpolated point acceptable
      if (2.0*p < 3.0*m*q - abs(tol*q)) & (p < abs(0.5*e*q))
         e = d;
         d = p/q;
      else
         d = m;
         e = m;
      end;
   end

   % Next point
   c = b;
   fc = fb;
   if abs(d) > tol
      b = b + d;
   else
      b = b - sign(b-a)*tol;
   end
   fb = F(b,varargin{:});
   k = k+1;
end
fprintf('Number of Iterations in fzerotx: %d\n',k)

Published with MATLAB® R2014a

function v = polyinterp(x,y,u)
%POLYINTERP  Polynomial interpolation.
%   v = POLYINTERP(x,y,u) computes v(j) = P(u(j)) where P is the
%   polynomial of degree d = length(x)-1 with P(x(i)) = y(i).

% Use Lagrangian representation.
% Evaluate at all elements of u simultaneously.

n = length(x);
v = zeros(size(u));
for k = 1:n
   w = ones(size(u));
   for j = [1:k-1 k+1:n]
      w = (u-x(j))./(x(k)-x(j)).*w;
   end
   v = v + w*y(k);
end

Published with MATLAB® R2014a