4

ニュートン法を 2 次元で実装する関数を作成しようとしていますが、これを行っている間に、関数の入力パラメーターが列ベクトル、ヤコビ行列の f(x) になるようにスクリプトを調整する必要があります。 of 、関数とそのヤコビ行列が別の .m ファイルにあるf(x)場合の初期推定と許容誤差。x0f(x)

私が書いたニュートン法を実装するスクリプトの例として、次のものがあります。

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

したがって、このスクリプトを使用して、関数とヤコビ行列を実際のスクリプトに実装しましたが、必要な入力パラメーターを使用して実際にスクリプトを作成する方法を見つけるのに苦労しています。

ありがとう！

4

1 に答える 1

7

よろしければ、コードを再構築して、より動的でユーザー フレンドリーに読みやすくしたいと思います。

いくつかの準備から始めましょう。スクリプトを真に動的にしたい場合は、Symbolic Math Toolbox を使用することをお勧めします。このように、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

xこれは 2D であるため、2D 関数が必要になります...と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

もう 1 つ... 方程式をベクトルまたは行列に配置し、 と のすべての値を各方程式に同時に代入するために使用できますsubs。xy

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]

行列に対しても同じことができますが、簡潔にするために、その方法は示しません。私はコードを延期します、そしてあなたはそれを見ることができます。

これで確立されたので、これを真に動的なものにするために、一度に 1 つずつコードに取り組みましょう。関数には、初期推定値、列ベクトルとしてx0の関数f(x)、2 x 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');

おそらく、Symbolic Math Toolbox を使用して計算を行っている場合、計算中の数値のデータ型はオブジェクトであることをお伝えしておく必要がありsymます。おそらくこれらを実数に変換し直したいので、 を使用doubleしてそれらをキャストバックできます。ただし、それらをsymフォーマットのままにしておくと、探しているものであれば、数値がきちんとした分数として表示されます。double小数点表現が必要な場合はにキャストします。

この関数を実行すると、探しているものが得られるはずです。このコードはまだテストしていませんが、うまくいくと確信しています。

ご不明な点がございましたら、お気軽にお問い合わせください。お役に立てれば。

乾杯！

于 2014-02-13T07:25:00.147 に答える