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

13 views (last 30 days)
Anthony Knighton
Anthony Knighton on 6 Mar 2021
Edited: Anthony Knighton on 7 Mar 2021
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
Anthony Knighton
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!

Translated by