C++ - 非線形方程式の解法(ニュートン法)!

Updated:


※この記事は12年以上前に投稿されたもので、情報が古い可能性があります。

少し前には、f(x)=0 の解を2分法により求める C++ アルゴリズムについて紹介しました。

今日は、方程式 f(x)=0 の解をニュートン法により求める C++ アルゴリズム についてです。

ニュートン法の概念・アルゴリズムは以下の通り。

NONLINEAR_EQ_NEWTON_01

以下、数式が多いので TEX で記載。

NONLINEAR_EQ_NEWTON_02

以下、C++ によるサンプルソースです。

0. 前提条件Permalink

  • Linux Mint 13 Maya (64bit) での作業を想定。
  • g++ (Ubuntu/Linaro 4.6.3-1ubuntu5) 4.6.3

1. C++ ソース作成Permalink

今回、検証に使用した方程式は、 x3x+1=0 で、グラフは以下のようになる。

NONLINEAR_EQ_NEWTON_03

そして、

  • C++ なのでオブジェクト指向な作りにしている。
  • 収束しない場合は最大50回ループするようにしている。

File: nonlinear_equation_newton.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
/*********************************************
 * 非線形方程式の解法 ( ニュートン法 )
 *********************************************/
#include <iostream>  // for cout
#include <math.h>    // for fabs()
#include <stdio.h>   // for printf()

// 方程式定義
#define f(x) (x * x * x - x + 1)
// f(x) の x における傾き ( f(x) を1回微分 )
#define g(x) (3 * x * x - 1)

using namespace std;

/*
 * 計算クラス
 */
class Calc
{
    // 各種定数
    static const double eps = 1e-08;  // 打ち切り精度
    static const int  limit = 50;     // 打ち切り回数

    // 各種変数
    double x, dx;  // x, dx 値
    int k;         // LOOP インデックス

    public:
        // 非線形方程式を解く(ニュートン法)
        void calcNonlinearEquation();
};

/*
 * 非線形方程式を解く(ニュートン法)
 */
void Calc::calcNonlinearEquation()
{
  // x 初期値設定
  x = -2.0;

  // 打ち切り回数 or 打ち切り誤差になるまで LOOP
  for (k = 1; k <= limit; k++) {
    dx = x;
    x = x - f(x) / g(x);
    if (fabs(x - dx) / fabs(dx) < eps) {
      printf("x=%f\n", x);
      break;
    }
  }

  // 収束しなかった場合
  if (k > limit)
    cout << "収束しない" << endl;
}

/*
 * メイン処理
 */
int main()
{
    try
    {
        // 計算クラスインスタンス化
        Calc objCalc;

        // 非線形方程式を解く(ニュートン法)
        objCalc.calcNonlinearEquation();
    }
    catch (...) {
        cout << "例外発生!" << endl;
        return -1;
    }

    // 正常終了
    return 0;
}

2. C++ ソースコンパイルPermalink

$ g++ nonlinear_equation_newton.cpp -o nonlinear_equation_newton

何も出力されなければ成功です。

3. 実行Permalink

$ ./nonlinear_equation_newton
x=-1.324718

値を方程式に代入して計算するとほぼ 0 になり、2分法で計算した結果とも等しいので、計算は正しいと考えてよいと思う。


数学はおもしろいけど、コンピュータで実証するというのもおもしろいです。

※ちなみに最近の当方の C++ アルゴリズムについての記事は、古い C によるアルゴリズムに関する書物を参考に C++ に移植した形態となっています。

以上。





 

Sponsored Link

 

Comments