使用MATLAB編寫實現牛頓法在兩個維度

Question

我想寫一個實現牛頓法在兩個維度和，而我已經這樣做了如下功能的功能，我到現在調整我的腳本，這樣的輸入參數我的功能必須在一個列向量被F（X），的f(x)雅可比矩陣，初始猜測x0和公差其中函數f(x)和其雅可比矩陣是在單獨的.m文件。

爲腳本我寫了一個實現牛頓方法的一個例子，我有：

n=0;   %initialize iteration counter  
eps=1;   %initialize error  
x=[1;1];  %set starting value 

%Computation loop  
while eps>1e-10&n<100 
    g=[x(1)^2+x(2)^3-1;x(1)^4-x(2)^4+x(1)*x(2)];   %g(x)  
    eps=abs(g(1))+abs(g(2));        %error  
    Jg=[2*x(1),3*x(2)^2;4*x(1)^3+x(2),-4*x(2)^3+x(1)]; %Jacobian  
    y=x-Jg\g;           %iterate  
    x=y;             %update x  
    n=n+1;            %counter+1  
end 

n,x,eps  %display end values 

這個腳本

所以，我已經實現的功能和雅可比矩陣爲實際的劇本，我很努力研究如何使用所需的輸入參數實際創建腳本。

謝謝！

Answer 1

如果你不介意的話，我想調整你的代碼，使其更具活力，更加用戶友好的閱讀。

讓我們先從一些準備工作。如果你想讓你的腳本真正動態，那麼我會建議你使用符號數學工具箱。這樣，您就可以使用MATLAB爲您處理函數的衍生物。您首先需要使用syms命令，然後使用所需的任何變量。這告訴MATLAB你現在將把這個變量當作「符號」（即不是常量）。讓我們先從一些基礎知識：

syms x; 
y = 2*x^2 + 6*x + 3; 
dy = diff(y); % Derivative with respect to x. Should give 4*x + 6; 
out = subs(y, 3); % The subs command will substitute all x's in y with the value 3 
        % This should give 2*(3^2) + 6*3 + 3 = 39 

因爲這是2D，我們將需要2D功能......所以，讓我們定義x和y作爲變量。該方法調用subs命令會稍有不同：

syms x, y; % Two variables now 
z = 2*x*y^2 + 6*y + x; 
dzx = diff(z, 'x'); % Differentiate with respect to x - Should give 2*y^2 + 1 
dzy = diff(z, 'y'); % Differentiate with respect to y - Should give 4*x*y + 6 
out = subs(z, {x, y}, [2, 3]); % For z, with variables x,y, substitute x = 2, y = 3 
           % Should give 56 

一件事...我們可以把公式爲向量或矩陣，並使用subs到x和y所有值同時代入每個方程。

syms x, y; 
z1 = 3*x + 6*y + 3; 
z2 = 3*y + 4*y + 4; 
f = [z1; z2]; 
out = subs(f, {x,y}, [2, 3]); % Produces a 2 x 1 vector with [27; 25] 

我們可以爲矩陣做同樣的事情，但爲了簡潔起見，我不會告訴你如何做到這一點。我會推遲到代碼，然後你可以看到它。

現在我們有建立，讓我們來解決你的代碼一塊在同一時間，真正使這一動態。您的功能需要初始猜測x0，作爲列向量的函數f(x)，雅可比矩陣作爲一個2×2矩陣和公差tol。

之前運行腳本，你將需要生成參數：

syms x y; % Make x,y symbolic 
f1 = x^2 + y^3 - 1; % Make your two equations (from your example) 
f2 = x^4 - y^4 + x*y; 
f = [f1; f2]; % f(x) vector 

% Jacobian matrix 
J = [diff(f1, 'x') diff(f1, 'y'); diff(f2, 'x') diff(f2, 'y')]; 

% Initial vector 
x0 = [1; 1]; 

% Tolerance: 
tol = 1e-10; 

現在，使你的腳本到一個函數：

% To run in MATLAB, do: 
% [n, xout, tol] = Jacobian2D(f, J, x0, tol); 
% disp('n = '); disp(n); disp('x = '); disp(xout); disp('tol = '); disp(tol); 

function [n, xout, tol] = Jacobian2D(f, J, x0, tol) 

% Just to be sure... 
syms x, y; 

% Initialize error 
ep = 1; % Note: eps is a reserved keyword in MATLAB 

% Initialize counter 
n = 0; 

% For the beginning of the loop 
% Must transpose into a row vector as this is required by subs 
xout = x0'; 

% Computation loop  
while ep > tol && n < 100 
    g = subs(f, {x,y}, xout); %g(x) 
    ep = abs(g(1)) + abs(g(2)); %error  
    Jg = subs(J, {x,y}, xout); %Jacobian 
    yout = xout - Jg\g;   %iterate  
    xout = yout;    %update x  
    n = n + 1;     %counter+1  
end 

% Transpose and convert back to number representation 
xout = double(xout'); 

我應該告訴你，當你」重新使用符號數學工具箱進行計算，計算它們時的數字數據類型爲sym對象。您可能想要將這些轉換回實數，因此您可以使用double將它們轉換回來。但是，如果您將它們留在sym格式中，它會將您的數字顯示爲整齊的分數，如果這正是您想要的。如果您想要小數點表示形式，請將其轉換爲double。

現在，當你運行這個函數時，它應該給你你正在尋找的東西。我沒有測試過這個代碼，但我很肯定這會起作用。

很高興回答您可能遇到的任何問題。希望這可以幫助。

乾杯！