連立非線形方程式のn次元ニュートンラフソン法による数値解法

以下の4元連立非線形方程式を解く。
f:id:seinzumtode:20200519082129p:plain

clear; close all;
syms f(x,y,u,v) g(x,y,u,v) h(x,y,u,v) k(x,y,u,v)
f(x,y,u,v) =  x^3-3*x*y^2 ...
              +3*((x^2-y^2)*u-2*x*y*v)...
              +3*((u^2-v^2)*x-2*u*v*y)...
              +u^3-3*u*v^2 ...
              +6*(x^2-y^2+2*(x*u-y*v)+u^2-v^2)...
              +21*(x+u)+32;
g(x,y,u,v) = 3*x^2*y-y^3 ...
              +3*(2*x*y*u+(x^2-y^2)*v) ...
              +3*(2*u*v*x+(u^2-v^2)*y) ...
              +(3*u^2*v-v^3) ...
              +6*(2*x*y+2*(y*u+x*v)+2*u*v) ...
              +21*(v+y);                   
h(x,y,u,v) = x^2-y^2-u^2+v^2-1;
m(x,y,u,v)= 2*x*y-2*u*v-1;
J = jacobian([f,g,h,m],[x,y,u,v]);
J1 = inv(J);

x=zeros(1,10);
y=zeros(1,10);
u=zeros(1,10);
v=zeros(1,10);

delta=1e-8;
res=1e4;
x(1)=-1;
y(1)=-2;
u(1)=-1;
v(1)=2;
k=1;
while res>delta^2
  next = [x(k);y(k);u(k);v(k)] ...
      - J1(x(k),y(k),u(k),v(k)) * ...
      [f(x(k),y(k),u(k),v(k));
      g(x(k),y(k),u(k),v(k));
      h(x(k),y(k),u(k),v(k));
      m(x(k),y(k),u(k),v(k))]   ;
  x(k+1) = next(1);
  y(k+1) = next(2);
  u(k+1) = next(3);
  v(k+1) = next(4);
  res = (f(x(k+1),y(k+1),u(k+1),v(k+1)))^2 ...
      + (g(x(k+1),y(k+1),u(k+1),v(k+1)))^2 ...
      + (h(x(k+1),y(k+1),u(k+1),v(k+1)))^2 ...
      + (m(x(k+1),y(k+1),u(k+1),v(k+1))^2);
  disp(eval(res));
  k=k+1;
end
fprintf('ans\n x = %f\n y = %f\n u = %f\n v = %f\n',x(k),y(k),u(k),v(k));

計算結果

ans
 x = -1.508467
 y = -0.189549
 u = -1.129368
 v = 0.189549