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

22 Ansichten (letzte 30 Tage)
Anthony Knighton am 6 Mär. 2021
Bearbeitet: Anthony Knighton am 7 Mär. 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 Kommentare4 ältere Kommentare anzeigen4 ältere Kommentare ausblenden
Walter Roberson am 6 Mär. 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 am 7 Mär. 2021
Bearbeitet: Anthony Knighton am 7 Mär. 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));

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Newton-Raphson Method finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by