绘制多个案例的 Newton-Raphson 解决方案的结果

Question

考虑以下问题：

我现在在这个问题的第三部分。我写了矢量循环方程（q=teta2、x=teta3 和 y=teta4）：

fval(1,1) = r2*cos(q)+r3*cos(x)-r4*cos(y)-r1;
fval(2,1) = r2*sin(q)+r3*sin(x)-r4*sin(y);

我有这2个函数，除了x和y之外的所有变量都给定了。我在 this video.

的帮助下找到了根源

现在我需要绘制 q 与 x 和 q 与 y 的图表，当 q 处于 [0,2pi] 时带有增量q为2.5度。我应该怎么做才能绘制图表？

以下是我目前的尝试：

function [fval,jac] = lorenzSystem(X)

%Define variables
x = X(1);
y = X(2);
q = pi/2;

r2 = 15
r3 = 50
r4 = 45
r1 = 40

%Define f(x)
fval(1,1)=r2*cos(q)+r3*cos(x)-r4*cos(y)-r1;
fval(2,1)=r2*sin(q)+r3*sin(x)-r4*sin(y);

%Define Jacobian
 jac = [-r3*sin(X(1)), r4*sin(X(2));
     r3*cos(X(1)), -r4*cos(X(2))];

%% Multivariate NR

%Initial conditions:
X0 = [0.5;1];
maxIter = 50;
tolX = 1e-6;

X = X0;
Xold = X0;
for i = 1:maxIter
    [f,j] = lorenzSystem(X);
    X = X - inv(j)*f;

    err(:,i) = abs(X-Xold);
    Xold = X;
    if (err(:,i)<tolX)
        break;
    end
end

Answer 1

请查看下面我的解决方案，并研究它与您自己的解决方案有何不同。

function [th2,th3,th4] = q65270276()
[th2,th3,th4] = lorenzSystem();

hF = figure(); hAx = axes(hF);
plot(hAx, deg2rad(th2), deg2rad(th3), deg2rad(th2), deg2rad(th4)); 
xlabel(hAx, '\theta_2')
xticks(hAx, 0:pi/3:2*pi);
xticklabels(hAx, {'[=10=]$','$\frac{\pi}{3}$','$\frac{2\pi}{3}$','$\pi$','$\frac{4\pi}{3}$','$\frac{5\pi}{3}$','\pi$'});
hAx.TickLabelInterpreter = 'latex';
yticks(hAx, 0:pi/6:pi);
yticklabels(hAx, {'[=10=]$','$\frac{\pi}{6}$','$\frac{\pi}{3}$','$\frac{\pi}{2}$','$\frac{2\pi}{3}$','$\frac{5\pi}{6}$','$\pi$'});
set(hAx, 'XLim', [0 2*pi], 'YLim', [0 pi], 'FontSize', 16);
grid(hAx, 'on');
legend(hAx, '\theta_3', '\theta_4')
end

function [th2,th3,th4] = lorenzSystem()
th2 = (0:2.5:360).';
[th3,th4] = deal(zeros(size(th2)));

% Define geometry:
r1 = 40;
r2 = 15;
r3 = 50;
r4 = 45;

% Define the residual:
res = @(q,X)[r2*cosd(q)+r3*cosd(X(1))-r4*cosd(X(2))-r1; ... Δx=0
             r2*sind(q)+r3*sind(X(1))-r4*sind(X(2))];     % Δy=0

% Define the Jacobian:
J = @(X)[-r3*sind(X(1)),  r4*sind(X(2));
          r3*cosd(X(1)), -r4*cosd(X(2))];

X0 = [acosd((45^2-25^2-50^2)/(-2*25*50)); 180-acosd((50^2-25^2-45^2)/(-2*25*45))]; % Accurate guess
maxIter = 500;
tolX = 1e-6;

for idx = 1:numel(th2)
  X = X0;
  Xold = X0;
  err = zeros(maxIter, 1); % Preallocation
  for it = 1:maxIter
    % Update the guess
    f = res( th2(idx), Xold );
    X = Xold - J(Xold) \ f;
%     X = X - pinv(J(X)) * res( q(idx), X ); % May help when J(X) is close to singular
    
    % Determine convergence
    err(it) = (X-Xold).' * (X-Xold);
    if err(it) < tolX
      break
    end
    % Update history
    Xold = X;    
  end
  
  % Unpack and store θ₃, θ₄
  th3(idx) = X(1);
  th4(idx) = X(2);
  
  % Update X0 for faster convergence of the next case:
  X0 = X;
end

end

几个注意事项：

所有计算均以度为单位。
我用的具体绘图代码没那么有趣，重要的是我提前定义了所有θ₂，然后遍历它们找到θ₃和θ₄（没有递归，就像在您自己的实现中所做的那样）。
第一种情况 (θ₂=0) 的初始猜测（实际上是解析解）可以通过使用以下定律手动解决问题（即“纸上”）找到cosines. The solver also works for other guesses, but you might need to increase maxIter. Also, for certain guesses (e.g. X(1)==X(2)), the Jacobian is ill-conditioned, in which case you can use pinv.
如果我的计算是正确的，这就是结果：

绘制多个案例的 Newton-Raphson 解决方案的结果

Plotting the results of a Newton-Raphson solution for multiple cases

matlab

plot

angle

newtons-method

numerical-methods