Newton's method iterations

I have a code written for Netwon's method below. However, not all initial starting guesses would converge. How would I find out which initial values for x would converge for the function in, say 20 or 30 or 50 iterations? For example, I'll start with x0 = 0. Does this converge? Then I'll try x0 = 0.1. Does this converge? The x0 = 0.2, then 0.3, etc. Which ones would converge?

I am trying to figure out how to summarize this information in a table. Below is the code that I have so far for Newton's Method.

%Newton's Method x0 = 7
f = @(x) 2*exp(-2*x) + 4*sin(x) - 2*cos(2*x);
fp = @(x) 4*(-exp(-2*x) + sin(2*x) + cos(x));
x0 = 7; 
N = 10; 
tol = 1E-6;
x(1) = x0;
n = 2; 
nfinal = N + 1; 
while (n <= N + 1)
  fe = f(x(n - 1));
  fpe = fp(x(n - 1));
  x(n) = x(n - 1) - fe/fpe;
  if (abs(fe) <= tol)
    nfinal = n; 
    break;
  end
  n = n + 1;
end
plot(0:nfinal - 1,x(1:nfinal),'o-')
title('Solution:')
xlabel('Iterations')
ylabel('X')
table([1:length(x)]',x')
x(n) = x(n - 1) - fe/fpe;
fprintf('%3d: %20g %20g\n', n, x(n), abs(fe));
if (abs(fe) <= tol)

Answer 1

[Alan Stevens]

I suggest you plot a graph of your function, then you can see where good initial estimates would be. For example

f = @(x) 2*exp(-2*x) + 4*sin(x) - 2*cos(2*x);
x = 0:0.1:20;
y = f(x);
plot(x,y),grid

[James Holden]

I do have a graph of the functions with its zeros mapped out already, over the interval that I am interested in.

f = @(x) 2*exp(-2*x) + 4*sin(x) - 2*cos(2*x);
fplot(f, [0 8])
hold on
y = 0;
fplot(y, [0 8], 'k')
plot(0, 0,'sr')
plot(2.768, 0,'sr')
plot(6.658, 0,'sr')
hold off

I wanted to find a way to test which starting values converge and which would not. For example, I know that 7 converges due to the code above. But what about 6.9 or 7.1 or what about 6.8 or 7.2? At some point, the initial values will be too far away and the method will not converge.

Is there a better way than just copy and pasting the code over and over again for each different value and looking at it manually? I was hoping to be able to create a table with intial values in increments of 0.1 in one column and display its value in another. Is this possible?

[Alan Stevens]

You shouldn't need to keep copying code. Just put everything in a loop along the lines of:

x0 = 0.1:0.1:8;

for j = 1:numel(x0)

    x(1) = x0(j);
    ... etc.
       
end

and keep a record inside the loop of which converge within N iterations and which don't.

[James Holden]

Okay, yeah. That seems to work. Now I have

x0 = 0.1:0.1:8;
for j = 1:numel(x0)
    f = @(x) 2*exp(-2*x) + 4*sin(x) - 2*cos(2*x);
    fp = @(x) 4*(-exp(-2*x) + sin(2*x) + cos(x));
    x(1) = x0(j);
    N = 50; 
    tol = 1e-6;
    n = 2; 
    nfinal = N + 1; 
    while (n <= N + 1)
      fe = f(x(n - 1));
      fpe = fp(x(n - 1));
      x(n) = x(n - 1) - fe/fpe;
      if (abs(fe) <= tol)
        nfinal = n; 
        break;
      end
      n = n + 1;
    end
end

Could i just table all of these with something like?

T = table(x0(:), n(:), fe(:))

Or how would I go about actually seeing the information and what not? The above is giving me the error "Error using table, All table variables must have the same number of rows."

[Alan Stevens]

More like this perhaps (note: I've set N to be 10 in order to get some fails; at N = 50 they all converge):

f = @(x) 2*exp(-2*x) + 4*sin(x) - 2*cos(2*x);
fp = @(x) 4*(-exp(-2*x) + sin(2*x) + cos(x));
x0 = 0.1:0.1:8;
cv = 0; fv = 0;
N = 10;
tol = 1e-6;

for j = 1:numel(x0)
    x = x0(j);
    n = 0;
    fe = 1;
    while (abs(fe) > tol) && (n <= N)
      n = n+1;
      fe = f(x);
      fpe = fp(x);
      x = x - fe/fpe;
    end
    
    if abs(fe)<=tol
       cv=cv+1;
       C(cv,:) = [x0(j), n];
    else
        fv = fv+1;
        F(fv,1) = x0(j);
    end
    
end

if cv>0
    disp('Converged')
    disp('[x0,  n]')
    disp(C)
end
disp('   ')
if fv>0
    disp(['Failed to converge in ' num2str(N) ' iterations'])
    disp('x0')
    disp(F)
end

[James Holden]

Thank you so much. That's almost exactly what I was looking for. Is there a way now to add all my fe values to the table as well? I tried

if abs(fe) <= tol
        cv = cv + 1;
        c(cv,:,:) = [x0(j), n, fe];
    else
        fv = fv + 1;
        F(fv, 1) = x0(j);
    end

but it's telling me that it couldn't perform the assignment because the size of the left is 1-by-2 where the size of the right is 1-by-3. Is the

c(cv,:,:)

part right to set up three collumns instead of just two?

[Alan Stevens]

Try this

f = @(x) 2*exp(-2*x) + 4*sin(x) - 2*cos(2*x);
fp = @(x) 4*(-exp(-2*x) + sin(2*x) + cos(x));
x0 = 0.1:0.1:8;
cv = 0; fv = 0;
N = 10;
tol = 1e-6;

for j = 1:numel(x0)
    x = x0(j);
    n = 0;
    fe = 1;
    while (abs(fe) > tol) && (n <= N)
      n = n+1;
      fe = f(x);
      fpe = fp(x);
      x = x - fe/fpe;
    end
    
    if abs(fe)<=tol
       cv=cv+1;
       feval1(cv,1) = fe;
       C(cv,:) = [x0(j), n];
    else
        fv = fv+1;
        feval2(fv,1) = fe;
        F(fv,1) = x0(j);
    end
    
end
format shorte
if cv>0
    disp('Converged')
    disp('[x0,  n,   fe]')
    disp([C, feval1])
end
disp('   ')
if fv>0
    disp(['Failed to converge in ' num2str(N) ' iterations'])
    disp('x0, fe')
    disp([F, feval2])
end

Answer 2

[Matt J]

https://www.mathworks.com/matlabcentral/answers/1465089-newton-s-method-for-root-approximation-using-multiple-initial-values?s_tid=srchtitle