MATLAB: Newton-Raphson method to determine roots of square root function
9 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/540421/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/540426/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/540431/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/540436/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/540426/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/540446/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/540451/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/540456/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/540461/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/540466/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/540471/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/540476/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/540481/image.png)
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/540486/image.png)
Velocity:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/540491/image.png)
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 comentarios
Walter Roberson
el 6 de Mzo. de 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
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)
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
el 7 de Mzo. de 2021
Editada: Anthony Knighton
el 7 de Mzo. de 2021
Respuestas (0)
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!