MATLAB: Newton-Raphson method to determine roots of square root function

Background: I am trying to implement the Newton-Raphson to determine the classical truning points of a particle in the potential . To simplify computation, I am normalizing and L as and , respectively. This way, I do not have to explicitly define and L in the code. Using this, the potential can now just be given by for the sake of simplicity. For an appropriate energy E, the particle oscillates between two turning points and . This occurs when , where is the velocity of the particle. From the conservation of energy, the velocity is given by . Thus, . The derivative of this expression is required for the Newton-Raphson method. I calculated the derivative analytically and found it to be . Is this calculation correct?

Potential and Velocity Graphs: First, graphed the in Mathematica to the determine appropriate range of E. Then, I defined an E and graphed the velocity to visually inspect where the roots are at that E. Initially, I chose E=0.5.

Potential:

Velocity:

Code:

The user inputs E and the brackets for the first and second turning points. m is just defined as 10. For the first run, I let E=.5, and I estimated the bracket for the first turning point as x1=.1 and x2=.2 and the bracket for the seond turning point as x1=.4 and x2=.5. I am having some problems with the code. I never seems to exit from the loop. It just continuously runs and does not actually calculate the turning points. Appologies for the longwinded question. I wanted to be thorough. Here is the MATLAB code.

clear all;
% Potential
% U(x)=(sin(2*pi*x)-(1/4)*sin(4*pi*x));
% Velocity
% v(x)=sqrt(2*(E-U(x))/m
% Energy
E=input('Enter Energy Value=>');
% Bracket for first turning point
x1=input('First Turning Point: Enter x1=>');
x2=input('First Turning Point: Enter x2=>');
m=10;
% define velocity function
v=@(x) sqrt((2*(E-sin(2*pi*x)+(1/4)*sin(4*pi*x)))/m);
% define derivative
dv=@(x) (pi*(cos(4*pi*x)-2*cos(2*pi*x)))/(sqrt(2)*sqrt(m)*sqrt((1/4)*sin(4*pi*x)-sin(2*pi*x)+E);
% control parameter
tol=1e-8;
% check bracket
for i=1:2
    v1=v(x1);
    v2=v(x2);
    if v1*v2>0
        error('bracket does not meet necessary condition');
    end
    % begin Newton-Raphson method
    xm=(x1+x2)/2;
    vm=v(xm);
    x=xm;
    vx=vm;
    n=0;
    while abs(vx)>tol
        dvx=dv(x);
        x=x-vx/dvx;
        vx=v(x);
        n=n+1;
    end
    xc(i)=x;
    % bracket for second turning point
    x1=input('Second Turning Point: Enter x1=>');
    x2=input('Second Turning Point: Enter x2=>');
end
fprintf('First Turning Point = %.8f\n',xc(1));
fprintf('Second Turning Point= %.8f\n',xc(2));

[Walter Roberson]

    while abs(vx)>tol
        dvx=dv(x);
        x=x-vx/dvx;
        vx=v(x);
        n=n+1;
    end

You do not have any fail-safe on that loop in case the calculation is diverging.

[Anthony Knighton]

How should I implement that?

[Walter Roberson]

MaxTries = 200;
converged = false;
for n = 1 : MaxTries
    if abs(vx)>tol
       converged = true;
       break;
    end
    dvx=dv(x);
    x=x-vx/dvx;
    vx=v(x);
    n=n+1;
end
if ~converged
    error('No convergence in %d iterations; time to get out the debugger!', MaxTries)
end

[Anthony Knighton]

Okay, well now the question becomes how to implement a root finding technique to find the root of a function that does not necessarily converge. Let's try it on an easier function, like for example. Obviously, the root is . How would you find this numerically using Newton-Raphson? Here is what I have tried so far. It just stays in the loop, so I assume that it does not converge:

clear all;
% Evaluate root of square root function using Newton-Raphson method.
% implement and test with trivial function (f=sqrt(x-a)), then try to test with
% complicated v(x) function.
% define function
a=2;
f=@(x) sqrt(x-a);
% define derivative
df=@(x)  1/(2*sqrt(x-a));
% begin Newton-Raphson method
tol=1e-8;
x0=2.1; % initial point
fx0=f(x0); % function value at initial point
x=x0;
fx=fx0;
n=0;
while abs(f(x))>tol
    dfx=df(x);
    x=x-fx/dfx;
    fx=f(x);
    n=n+1;
end
xc=x
fprintf('root =%.8f\n',xc)

[Walter Roberson]

You have problems ;-)

% Evaluate root of square root function using Newton-Raphson method.
% implement and test with trivial function (f=sqrt(x-a)), then try to test with
% complicated v(x) function.
% define function
a=2;
f=@(x) sqrt(x-a);
% define derivative
df=@(x)  1/(2*sqrt(x-a));
% begin Newton-Raphson method
tol=1e-8;
x0=2.1; % initial point
fx0=f(x0); % function value at initial point
x=x0;
fx=fx0;
MaxTries = 200;
converged = false;
fx_record = zeros(MaxTries,1);
for n = 1 : MaxTries
    fx_record(n) = fx;
    if abs(fx)<tol
       converged = true;
       break;
    end
    dfx=df(x);
    x=x-fx/dfx;
    fx=f(x);
end

if ~converged
    fprintf('No convergence in %d iterations; time to get out the debugger!', MaxTries)
else
    xc=x;
    fprintf('root =%.8f\n',xc);
    fprintf('iterations = %d\n', n);
end

No convergence in 200 iterations; time to get out the debugger!

fx_record = fx_record(1:n);    %trim unused entries
fx_is_real_idx = find(imag(fx_record) == 0);
real_record = fx_record(fx_is_real_idx);
[~,idx] = min(abs(real_record));
closest_to_zero = real_record(idx)

closest_to_zero = 0.3162

subplot(2,1,1); plot(1:n, real(fx_record), 'k', idx, closest_to_zero, 'b*'); title('real(fx)')

subplot(2,1,2); plot(1:n, imag(fx_record), 'r'); title('imag(fx)')

[Anthony Knighton]

I figured out a simple way to determine turning points just using potential. The turning points occur where , so you can just solve using bisection or Newton-Raphson method. As you can see, that also solves . It keeps you from having to deal with the square root function. I don't know why I didn't realize this immediately. Here is the code if anyone is interested:

clear all;
% Potential
% U(x)=(sin(2*pi*x)-(1/4)*sin(4*pi*x));
% Velocity
% v(x)=sqrt(2*(E-U(x))/m
% Energy
E=input('Enter Energy Value=>');
% Bracket for first turning point
x1=input('First Turning Point: Enter x1=>');
x2=input('First Turning Point: Enter x2=>');
% define velocity function
U=@(x) sin(2*pi*x)-(1/4)*sin(4*pi*x)-E;
% define derivative
dU=@(x) 2*pi*cos(2*pi*x)-pi*cos(4*pi*x);
% control parameter
tol=1e-8;
% check bracket
for i=1:2
    U1=U(x1);
    U2=U(x2);
    if U1*U2>0
        error('bracket does not meet necessary condition');
    end
    % begin bisection method
    n=0;
    xm=(x1+x2)/2;
    Um=U(xm);
    while n<10;
        if U1*Um<0;
            x2=xm;
            U2=Um;
        else
            U1=xm;
            U2=Um;
        end
        xm=(x1+x2)/2;
        Um=U(xm);
        n=n+1;
    end
    fprintf('Root by bisection method     = %.8f  (iteration= %i)\n',xm,n)
    % begin Newton-Raphson method
    Um=U(xm);
    x=xm;
    Ux=Um;
    n=0;
    while abs(Ux)>tol
        dUx=dU(x);
        x=x-Ux/dUx;
        Ux=U(x);
        n=n+1;
    end
    fprintf('Root by Newton-Raphson method = %.8f  (iteration= %i)\n',x,n);
    xc(i)=x;
    % bracket for second turning point
    x1=input('Second Turning Point: Enter x1=>');
    x2=input('Second Turning Point: Enter x2=>');
end
fprintf('First Turning Point = %.8f\n',xc(1));
fprintf('Second Turning Point= %.8f\n',xc(2));