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

Question

Недавно преподавателю был дан фрагмент кода Matlab для решения одновременных уравнений с использованием метода Ньютона-Рафсона с матрицей якобиан (I В комментариях он также остался. Однако, хотя он предоставил мне базовый код, я не могу заставить его работать, как бы я ни старался. Я потратил много часов, пытаясь представить функцию «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).

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

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.

---

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

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