以下の内容はhttps://seinzumtode.hatenadiary.jp/entry/2020/05/19/065758より取得しました。


2次元ニュートンラフソン法による3次方程式の複素解の数値解法

 x^3+6x^2+21x+32=0
カルダノの方法によって、解は
 x_1=-\sqrt[3]{9}+\sqrt[3]{3}-2=-2.637834253
 x_2= -\sqrt[3]{9}e^{\frac{2}{3}\pi i}+\sqrt[3]{3}e^{\frac{4}{3}\pi i}-2=1.6810829-3.0504302i
 x_3= -\sqrt[3]{9}e^{\frac{4}{3}\pi i}+\sqrt[3]{3}e^{\frac{2}{3}\pi i}-2=-1.6810829‐ +3.0504302i

clear; close all;

f = @(x,y) x^3-3*x*y^2+6*x^2-6*y^2+21*x+32;
g = @(x,y) 3*x^2*y-y^3+12*x*y+21*y;
dfdx = @(x,y) 3*x^2-3*y^2+12*x+21;
dfdy = @(x,y) -6*x*y-12*y;
dgdx = @(x,y) 6*x*y+12*y;
dgdy = @(x,y) 3*x^2-3*y^2+12*x+21;

x=zeros(1,100);
y=zeros(1,100);
x(1)=-1;
y(1)=-2;

delta=1e-10;
res=1e4;
k=2;
while res>delta^2
  J = dfdx(x(k-1),y(k-1))*dgdy(x(k-1),y(k-1))-...
      dfdy(x(k-1),y(k-1))*dgdx(x(k-1),y(k-1));
  x(k) = x(k-1)-1/J*...
      (f(x(k-1),y(k-1))*dgdy(x(k-1),y(k-1))-...
       g(x(k-1),y(k-1))*dfdy(x(k-1),y(k-1)))  ;
  y(k) = y(k-1)-1/J*...
      (g(x(k-1),y(k-1))*dfdx(x(k-1),y(k-1))-...
       f(x(k-1),y(k-1))*dgdx(x(k-1),y(k-1)))  ;
  res=f(x(k),y(k))^2+g(x(k),y(k))^2;
  %disp(res);
  k = k+1;
end
fprintf('ans\n z= (%d) + (%d) i\n',x(k-1),y(k-1));

計算結果

ans
 z= (-1.681083e+00) + (-3.050430e+00) i

複素解の一つ x_2が計算できている。
もう一つの複素解 x_3 x_2の共役となる。




以上の内容はhttps://seinzumtode.hatenadiary.jp/entry/2020/05/19/065758より取得しました。
このページはhttp://font.textar.tv/のウェブフォントを使用してます

不具合報告/要望等はこちらへお願いします。
モバイルやる夫Viewer Ver0.14