Why doesn't my code work for Newton's Method

Find a minimum using the Newton's method with .

f = @(x) cos(x) ./ x;
x0 = 3;
N = 10;
%below is my code, f_p is f prime and f_dp is f double prime
function min=Newtons_Method(f,x0,N)
dx=1e-4;
x1=x0+dx;
f_p=@(x) (f(x1)-f(x0))/dx;
f_dp=@(x) (f_p(x1)-f_p(x0))/dx;
for i=1:N
    x2=x0-(f_dp(x0)/f_p(x0));
    x0=x2;
    x1=x0+dx;
    f=f(x0);
    f_p=@(x) (f(x1)-f(x0))/dx;
    f_dp=@(x) (f_p(x1)-f_p(x0))/dx;
end
min=x0;
end

[Walter Roberson]

Answer 1

[Walter Roberson]

f = @(x) cos(x) ./ x;
x0 = 3;
N = 10;
%below is my code, f_p is f prime and f_dp is f double prime

min_found_at = Newtons_Method(f, x0, N)

Array indices must be positive integers or logical values.

Error in solution>@(x)(f(x1)-f(x0))/dx (line 18)
f_p=@(x) (f(x1)-f(x0))/dx;

Error in solution>@(x)(f_p(x1)-f_p(x0))/dx (line 19)
f_dp=@(x) (f_p(x1)-f_p(x0))/dx;

Error in solution>Newtons_Method (line 14)
x2=x0-(f_dp(x0)/f_p(x0));

function min=Newtons_Method(f,x0,N)
dx=1e-4;
x1=x0+dx;
f_p=@(x) (f(x1)-f(x0))/dx;
f_dp=@(x) (f_p(x1)-f_p(x0))/dx;
for i=1:N
    x2=x0-(f_dp(x0)/f_p(x0));
    x0=x2;
    x1=x0+dx;
    f=f(x0);
    f_p=@(x) (f(x1)-f(x0))/dx;
    f_dp=@(x) (f_p(x1)-f_p(x0))/dx;
end
min=x0;
end

First consider

f_p=@(x) (f(x1)-f(x0))/dx;

That function handle ignores x and calculates f (passed in) at the original location x0 passed in, and at x0 (passed in) plus dx . f_p is not a function of x0 or x1 and it ignores x, so if x0 or x1 or x were to change, it would ignore them and calculate at the fixed locations.

f_dp=@(x) (f_p(x1)-f_p(x0))/dx;

since f_p ignores its inputs, f_p(x1) - f_p(x0) is going to be "some value" minus the same value. So f_dp is going to return 0.

for i=1:N
    x2=x0-(f_dp(x0)/f_p(x0));

first iteration, that is going to be 0 divided by some fixed value, and subtract the resulting 0 from x0, so x2 is going to be the same as x0

    x0=x2;
    x1=x0+dx;

Yes, but that does not change x1 or x0 that are already written into the function handles.

    f=f(x0);

First iteration, you evaluate the passed-in f at the (updated) x0 getting a numeric result. And you store that numeric result into f, overwriting the function handle you passed in.

    f_p=@(x) (f(x1)-f(x0))/dx;

This overwrites f_p with a new function handle. Execution of that function handle would try to index the scalar numeric value now stored in f at the fixed location x0 and x1. That going to fail because x0 and x1 are likely not valid indices into the scalar numeric value now stored in f.

If somehow it did succeed, the result would be independent of x and would use whatever x0 and x1 are present the workspace at the time the function handle is defined, ignoring any change to x0 and x1.

    f_dp=@(x) (f_p(x1)-f_p(x0))/dx;

You overwrite function handle f_p with new function handle f_p on the previous line, so this line is syntactically valid... but suffers from the problem that f_p ignores its input so the function handle is always going to calculate 0 (unless the original f generated inf or -inf or nan or empty)

[Walter Roberson]

As you are doing a numeric approximation, you might want to consider something like

f_p = @(x) (f(x+dx/2) - f(x-dx/2)) ./ dx

and that function handle would not need to be updated.

[Ryley]

Thank you, this cleared up a lot for me.