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

26 views (last 30 days)
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));
  6 Comments
Walter Roberson
Walter Roberson on 6 Mar 2021
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
Anthony Knighton on 7 Mar 2021
Edited: Anthony Knighton on 7 Mar 2021
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));

Sign in to comment.

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!