A..3 complex-newton.C

/*
 * complex-newton.C -- Newton 法で方程式 f(x)=0 を解く
 *  コンパイル: g++ -o complex-newton complex-newton.C
 *  実行: ./complex-newton
 *
 *  2016/3/22修正
 */

#include <iostream>
#include <complex>
using namespace std;

int main(void)
{
    int i, maxitr = 100;
    complex<double> f(complex<double>), dfdz(complex<double>), x, dx;
    double eps;

    cout << " 初期値 x0, 許容精度ε=";
    cin >> x >> eps;

    cout.precision(16);
    for (i = 0; i < maxitr; i++) {
        dx = - f(x) / dfdz(x);
        x += dx;
        cout << "f(" << x << ")=" << f(x) << endl;
        if (abs(dx) <= eps) break;
    }
    return 0;
}

/*
 * この関数の例は杉原・室田『数値計算法の数理』岩波書店, p.67 による
 *   f(z) = z^3-2z+2
 *   1実根(≒-1.769292354238631),2虚根(≒0.8846461771193157±0.5897428050222054)
 *   を持つが、原点の近くに初期値を取ると、0 と 1 の間を振動する。
 */

complex<double> f(complex<double> z)
{
  /* return z * z * z - 2 * z + 2; */
  return z * (z * z - 2.0) + 2.0;
}

/* 関数 f の導関数 (df/dx のつもりで名前をつけた) */
complex<double> dfdz(complex<double> z)
{
  return 3.0 * z * z - complex<double>(2,0);
}

コンパイル＆実行例

% g++ complex-newton.C
% ./a.out
 初期値 x0, 許容精度ε=0.3 1e-14
f((1.12485549132948,0))=(1.173568531458942,0)
f((0.4713843772268934,0))=(1.161974377253067,0)
f((1.342827923349606,0))=(1.735713781960667,0)
f((0.8337553363070073,0))=(0.9120726492429374,0)
f((-9.840767203756229,0))=(-931.3052418605159,0)
f((-6.612920012816673,0))=(-273.9618545552224,0)
f((-4.4923431047393,0))=(-79.67594843313302,0)
f((-3.131371682066463,0))=(-22.44188600195348,0)
f((-2.312816652470607,0))=(-5.745902514520807,0)
f((-1.903778836863019,0))=(-1.092448577226583,0)
f((-1.780659980946024,0))=(-0.0847076152532642,0)
f((-1.769384049454865,0))=(-0.0006777810558658004,0)
f((-1.769292360276141,0))=(-4.462436198338082e-08,0)
f((-1.769292354238631,0))=(2.220446049250313e-16,0)
f((-1.769292354238631,0))=(2.220446049250313e-16,0)
mk@katsuradanoMacBook-Pro-555 progs-algebraic-equation % ./a.out
 初期値 x0, 許容精度ε=0.1 1e-14
f((1.014213197969543,0))=(1.014822114238246,0)
f((0.07965576631987636,0))=(1.841193886472037,0)
f((1.009098740372765,0))=(1.009347854859992,0)
f((0.05222652653371329,0))=(1.895689400532465,0)
f((1.003965184727484,0))=(1.004012415140623,0)
f((0.02332943565497392,0))=(1.953353826028611,0)
f((1.000804353182403,0))=(1.000806294654933,0)
f((0.004806794546775128,0))=(1.990386521968754,0)
f((1.000034548045781,0))=(1.000034551626525,0)
f((0.0002072524744252124,0))=(1.999585495060052,0)
f((1.000000064421484,0))=(1.000000064421497,0)
f((3.865287804272199e-07,0))=(1.999999226942439,0)
f((1.000000000000224,0))=(1.000000000000224,0)
f((1.34559030584569e-12,0))=(1.999999999997309,0)
f((1,0))=(1,0)
f((0,0))=(2,0)
f((1,0))=(1,0)
f((0,0))=(2,0)

中略

f((0,0))=(2,0)
f((1,0))=(1,0)
f((0,0))=(2,0)
%

初期値 $0.3$ で始めると実根に収束するが、 初期値 $0.1$ で始めると、途中から $(0,0)$ と $(1,0)$ の間を振動する。 確かに本に書いてある通り。