Как использовать метод Ньютона-Рафсона в Matlab, чтобы найти корень уравнения?

Question

Я новый пользователь MATLAB. Я хочу найти значение, которое составляет f(x) = 0, используя метод Ньютона-Рафсона. Я пытался написать код, но кажется, что трудно реализовать метод Ньютона-Рафсона. Это то, что я до сих пор:

function x = newton(x0, tolerance)
    tolerance = 1.e-10;
    format short e;
    Params = load('saved_data.mat');
    theta = pi/2;
    zeta = cos(theta);
    I = eye(Params.n,Params.n);
    Q = zeta*I-Params.p*Params.p';

    % T is a matrix(5,5)
    Mroot = Params.M.^(1/2);  %optimization
    T = Mroot*Q*Mroot;

    % Find the eigenvalues
    E = real(eig(T));

    % Find the negative eigenvalues
    % Find the smallest negative eigenvalue
    gamma = min(E);  

    % Now solve for lambda
    M_inv = inv(Params.M);  %optimization
    zm = Params.zm;

    x = x0;
    err = (x - xPrev)/x;

    while abs(err) > tolerance
        xPrev = x;
        x = xPrev - f(xPrev)./dfdx(xPrev);

        % stop criterion: (f(x) - 0) < tolerance
        err = f(x);
   end 

   % stop criterion: change of x < tolerance % err = x - xPrev;

end

Вышеуказанная функция используется так:

% Calculate the functions
Winv = inv(M_inv+x.*Q);

f = @(x)( zm'*M_inv*Winv*M_inv*zm);

dfdx = @(x)(-zm'*M_inv*Winv*Q*M_inv*zm);

x0 = (-1/gamma)/2;

xRoot = newton(x0,1e-10);

Answer 1

Вопрос не особо понятен. Тем не менее, вам нужно реализовать рут найти себя? Если нет, то просто используйте встроенную функцию Matlab fzero (не основанную на Ньютоне-Рафсоне).

Если вам нужна собственная реализация метода Ньютона-Рафсона, тогда я предлагаю использовать один из ответов на метод Ньютона-Рафсона в Matlab? в качестве отправной точки.

Редактировать: Следующее не отвечает на ваш вопрос, но просто примечание о стиле кодирования.

Полезно разбить вашу программу на куски многократного использования. В этом случае ваш корень должен быть отделен от конструкции вашей функции. Я рекомендую записать ваш метод Ньютона-Рафсона в отдельный файл и вызвать его из сценария, в котором вы определяете свою функцию и ее производную. Ваш источник будет выглядеть примерно так:

% Define the function (and its derivative) to perform root finding on:
Params = load('saved_data.mat');
theta = pi/2;
zeta  = cos(theta);
I = eye(Params.n,Params.n);
Q = zeta*I-Params.p*Params.p';

Mroot = Params.M.^(1/2);
T = Mroot*Q*Mroot; %T is a matrix(5,5)

E = real(eig(T)); % Find the eigen-values

gamma = min(E);   % Find the smallest negative eigen value

% Now solve for lambda (what is lambda?)
M_inv = inv(Params.M);
zm = Params.zm;

Winv = inv(M_inv+x.*Q);

f = @(x)( zm'*M_inv*Winv*M_inv*zm);
dfdx = @(x)(-zm'*M_inv*Winv*Q*M_inv*zm);

x0 = (-1./gamma)/2.;

xRoot = newton(f, dfdx, x0, 1e-10);

В newton.m у вас будет реализация метода Ньютона-Рафсона, который принимает в качестве аргументов определяемые вами функции (f и dfdx). Используя ваш код, указанный в вопросе, это будет выглядеть примерно так:

function root = newton(f, df, x0, tol)

    root = x0; % Initial guess for the root

    MAXIT = 20; % Maximum number of iterations

    for j = 1:MAXIT;

        dx = f(root) / df(root);
        root = root - dx

        % Stop criterion:
        if abs(dx) < tolerance
            return
        end

    end

    % Raise error if maximum number of iterations reached.
    error('newton: maximum number of allowed iterations exceeded.')

end 

Обратите внимание, что я избежал использования бесконечного цикла.