Пытаясь решить уравнения одновременности в Matlab, не может понять, как форматировать функции

Question

Недавно один лектор дал мне фрагмент кода Matlab за способ решения одновременных уравнений с использованием метода Ньютона-Рафсона с помощью матрицы Якоби (я также оставил в его комментариях). Тем не менее, хотя он и предоставил мне основной код, я не могу заставить его работать, как бы я ни старался. Я потратил много часов, пытаясь представить функцию "func", но безрезультатно, часто получая сообщение о том, что не хватает входных данных. Буду признателен за любую помощь, особенно в том, как написать функцию 'func'.

  function root = newtonRaphson2(func,x,tol)

% Newton-Raphson method of finding a root of simultaneous    
% equations fi(x1,x2,...,xn) = 0, i = 1,2,...,n.    
% USAGE: root = newtonRaphson2(func,x,tol)    
% INPUT:    
% func = handle of function that returns[f1,f2,...,fn].
% x = starting solution vector [x1,x2,...,xn].
% tol = error tolerance (default is 1.0e4*eps).
% OUTPUT:
% root = solution vector.

if size(x,1) == 1; x = x'; end % x must be column vector    
for i = 1:30    
   [jac,f0] = jacobian(func,x);    
   if sqrt(dot(f0,f0)/length(x)) < tol   
      root = x; return
   end   

   dx = jac\(-f0);   
   x = x + dx;    
   if sqrt(dot(dx,dx)/length(x)) < tol    
      root = x; return    
   end
end

error('Too many iterations')

function [jac,f0] = jacobian(func,x)

% Returns the Jacobian matrix and f(x).

h = 1.0e-4;    
n = length(x);    
jac = zeros(n);    
f0 = feval(func,x);    
for i =1:n    
   temp = x(i);    
   x(i) = temp + h;    
   f1 = feval(func,x);    
   x(i) = temp;    
   jac(:,i) = (f1 - f0)/h;   
end

Решаемые одновременные уравнения:

sin(x)+y^2+ln(z)-7=0
3x+2^y-z^3+1=0
x+y+Z-=0

с отправной точкой (1,1,1).

Тем не менее, они являются произвольными и могут быть заменены на что угодно, в основном мне просто нужно знать общий формат. Большое спасибо, я знаю, что это может быть очень простой задачей, но я только недавно начал учить себя матлабу.

Answer 1

Вам нужно создать новый файл с именем myfunc.m (или любое другое имя), которое принимает один входной параметр - вектор столбца x - и возвращает один выходной вектор - вектор столбца y такой, что y = f(x),

Например,

function y = myfunc(x)

    y = zeros(3, 1);
    y(1) = sin(x(1)) + x(2)^2 + log(x(3)) - 7;
    y(2) = 3*x(1) + 2^x(2) - x(3)^3 + 1;
    y(3) = x(1) + x(2) + x(3);

end

Затем вы можете обратиться к этой функции как @myfunc как в

>> newtonRaphson2(@myfunc, [1;1;1], 1e-6);

Причина специальной записи в том, что Matlab позволяет вам вызывать функцию без параметров, опуская символы () которые следуют за этим. Так, например, Matlab интерпретирует myfunc когда вы вызываете функцию без аргументов (поэтому она пытается заменить ее своим результатом), тогда как @myfunc относится к самой функции, а не к ее результату.

В качестве альтернативы вы можете написать функцию напрямую, используя @ обозначение, как в

>> newtonRaphson2(@(x) exp(x) - 3*x, 2, 1e-2)
ans =
     1.5315
>> newtonRaphson2(@(x) exp(x) - 3*x, 1, 1e-2)
ans =
     0.6190

которые являются двумя корнями уравнения exp(x) - 3 * x = 0,

---

Отредактируйте - кроме этого, у вашего профессора ужасный стиль кодирования (если код в вашем вопросе является прямой копией того, что он вам дал, и вы не искали его по пути). Было бы лучше написать код, подобный этому, с отступом, поясняющим, какова структура кода.

function root = newtonRaphson2(func, x, tol)
% Newton-Raphson method of finding a root of simultaneous
% equations fi(x1,x2,...,xn) = 0, i = 1,2,...,n.
%
% USAGE: root = newtonRaphson2(func,x,tol)
%
% INPUT:
%   func = handle of function that returns[f1,f2,...,fn].
%   x    = starting solution vector [x1,x2,...,xn].
%   tol  = error tolerance (default is 1.0e4*eps).
%
% OUTPUT:
%   root = solution vector.

    if size(x, 1) == 1; % x must be column vector
        x = x';
    end 

    for i = 1:30
        [jac, f0] = jacobian(func, x);

        if sqrt(dot(f0, f0) / length(x)) < tol
            root = x; return
        end

        dx = jac \ (-f0);
        x  = x + dx;

        if sqrt(dot(dx, dx) / length(x)) < tol
            root = x; return
        end
    end

    error('Too many iterations')

end

function [jac, f0] = jacobian(func,x)
% Returns the Jacobian matrix and f(x).

    h = 1.0e-4;
    n = length(x);
    jac = zeros(n);
    f0 = feval(func,x);

    for i = 1:n
        temp = x(i);
        x(i) = temp + h;
        f1 = feval(func,x);
        x(i) = temp;
        jac(:,i) = (f1 - f0)/h;
    end

end