如果您不介意,我想重组您的代码,使其更加动态且更易于阅读。
让我们从一些预备知识开始。如果您想让脚本真正动态化,那么我建议您使用符号数学工具箱。这样,您就可以使用 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
and 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
and 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 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');
我可能应该告诉您,当您使用符号数学工具箱进行计算时,计算时数字的数据类型是sym
目的。您可能想将它们转换回实数,因此您可以使用double
把他们赶回去。但是,如果您将它们留在sym
格式,如果您正在寻找的话,它会将您的数字显示为整齐的分数。投射到double
如果你想要小数点表示形式。
现在,当您运行此函数时,它应该会为您提供所需的内容。我还没有测试过这段代码,但我很确定这会起作用。
很高兴回答您可能有的更多问题。希望这可以帮助。
Cheers!