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

Question

Anthony Knighton 2021-3-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/764386-matlab-newton-raphson-method-to-determine-roots-of-square-root-function

编辑： Anthony Knighton 2021-3-7

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 个评论
显示 4更早的评论隐藏 4更早的评论

Walter Roberson 2021-3-6

在 MATLAB Online 中打开

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 2021-3-7

编辑：Anthony Knighton 2021-3-7

在 MATLAB Online 中打开

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));

请先登录，再进行评论。

请先登录，再回答此问题。

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

6 个评论
显示 4更早的评论隐藏 4更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

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

6 个评论 显示 4更早的评论隐藏 4更早的评论

回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

6 个评论
显示 4更早的评论隐藏 4更早的评论