プログラミング演習1


― 課題4 数値解析 ―

2.2 ニュートン法

ニュートン法は以下の手順で方程式 f (x)=0 の解を求める。

  1. 初期値x0を与える。
  2. その点での関数値f (x0)と 傾きf '(x0)から x 切片x1を推定する。
  3. x1を新たな初期値x0としてこの操作を反復する。

(x0, f (x0))を通る接線の式
　　y - f (x0) = f '(x0) (x - x0)
においてx 切片では y = 0, x = x1 となるので
　　x1 = x0 - f (x0) / f '(x0)
f (x) = x2 - 2 のとき、f '(x) = 2x だから
　　x1 = x0 - (x02 - 2) / 2x0 = (x0 + 2 / x0) / 2

真野芳久「Pascalプログラミングの基礎」サイエンス社 p.41 のプログラムと 川上一郎「数値計算」岩波書店 p.180 付録3.1のプログラムを合体してCに移植し、 手直ししたもの(数値解析の講義用Webページで公開している Newton.c を一部手直し) を以下に示す。
# ただし、練習のためにバグが入れてある。
下線部では scanf関数の代わりにgets関数とsscanf関数 を使って入力する方法を用いた。
また、リダイレクト出力時にも入力促進メッセージ等がコマンドプロンプトに 出力されるようにfprintf関数を用いて 標準エラー出力(stderr)とした。 

 1: /*
 2:  * Newton.c: ニュートン法
 3:  */
 4: #include <stdio.h> // printf, gets, sscanf
 5: #include <math.h>  // fabs
 6:
 7: main()
 8: {
 9:     int n = 0;
10:     double a, x, x0, err, eps = 1.0 e-10;
11:     char s[128];
12:
13:     fprintf(stderr, " a = "); gets(s); sscanf(s, "%lf", &a);
14:     while (a <= 0.0) {
15:         fprinrf(stderr, "'a' には正の数を入れてください。\n");
16:         fprintf(stderr, " a = "); gets(s); sscanf(s, "%1f", &a)
17:     }
18:
19: 　　x = (a + 1.0) / 2.0;
20:     printf("# n, x, err\n");
21:     printf("%4d, % .15e\n", n, x);
22:     do {
23:         n ++;
24:         x0 = x;
25:         x = (x0 + a/x0) / 2;
26:         err = fabs(x - x0);
27:         printf("%4d, % .l5e, % .l5e\n", n, x, err)
28:     } while (err >= eps);
29:     printf("\n# sqrt(%e) =  %.15e\n"，a, x);
30:     retrun (O);
31: }

主要部分のアルゴリズムをPADで描くと

" a = 2" を入力したときの実行結果とそのグラフは

2.3 二元連立非線形方程式（ニュートン法）

n 元連立非線形方程式は一般に
　　f1(x1, x2, ..., xn) = 0,
　　f2(x1, x2, ..., xn) = 0,
　　　...　...　...　...　...
　　fn(x1, x2, ..., xn) = 0
と書くことができる。これをベクトル表記すると
　　f (x) = 0
となる。

このとき、ニュートン法の計算式は
　　xk+1 = xk + u (xk)
　　　　　 = xk - J (xk)-1 f (xk)
となる。ここに、J (xk) はヤコビ行列と呼ばれ、
　　
であり、また、 u (x) = - J (x)-1 f (x) である。

2変数の場合は、逆行列が簡単に求められるので
　　
ここに
　　

「数値解析」の講義用Webページで公開している 「連立非線形方程式」を解くプログラム(Newton2.c)を 一部手直ししたものを以下に示す。 このプログラムは
　　x 2 + y 2 = 1, 　　　 y = x 3
を解くためのものである。仕様は以下のとおり。

  1. 入力は初期値(x0, y0)のみ。許容誤差は eps = 10-10。
  2. 出力は途中経過と最終結果。CSV形式に対応。
  3. 方程式の左辺の関数2つはおのおの独立に定義。（結果は戻り値で返す。）
  4. ヤコビ行列の4つの要素もおのおの独立に定義。（結果は戻り値で返す。）
  5. 補正ベクトルu (x)計算ルーチンの結果は2つの引数で返す。
  6. 2変数ニュートン法ルーチンの計算結果は2つの引数で返す。 
     1: /*
 2:  *  Newton2.c: 2変数のニュートン法
 3:  *  x2+y2=1, y=x3 -> (x,y)=(0.826,0.5636),(-0.826,-0.5636)
 4:  */
 5: #include <stdio.h> // printf, fprintf, gets, sscanf
 6: #include <math.h>  // hypot(hypotenuse): 直角三角形の斜辺を求める関数
 7:
 8: double f1(double x, double y)
 9: {
10:     return x * x + y * y - 1; //  f1(x,y) = x2 + y2 - 1
11: }
12: double f2(double x, double y)
13: {
14:     return x * x * x - y;     //  f2(x,y) = x3 - y
15: }
16: 
17: double J11(double x, double y)
18: {
19:     return 2 * x;     //  df1(x,y)/dx = 2x
20: }
21: double J12(double x, double y)
22: {
23:     return 2 * y;     //  df1(x,y)/dy = 2y
24: }
25: double J21(double x, double y)
26: {
27:     return 3 * x * x; //  df2(x,y)/dx = 3x2
28: }
29: double J22(double x, double y)
30: {
31:     return -1;        //  df2(x,y)/dy = -1
32: }
33: 
34: int ucalc(double x, double y, double *ux, double *uy)
35: /*
36:  *  u = -J-1(x)*f(x)
37:  */
38: {
39:     double det;
40: 
41:     det = J11(x, y) * J22(x, y) - J12(x, y) * J21(x, y);
42:     *ux = -1 / det * (J22(x, y) * f1(x, y) - J12(x, y) * f2(x, y));
43:     *uy =  1 / det * (J21(x, y) * f1(x, y) - J11(x, y) * f2(x, y));
44:     return 0;
45: }
46: 
47: int Newton2(double *x, double *y, double eps)
48: {
49:     int n=0;
50:     double x0, y0, ux, uy, err;
51: 
52:     printf("# n, x, y, err\n");
53:     printf("%4d, % .15e, % .15e\n", n, *x, *y);
54:     do {
55:         n++;
56:         x0 = *x; y0 = *y;
57:         ucalc(*x, *y, &ux, &uy);
58:         *x += ux; *y += uy;        // xk+1 = xk + u(xk)
59:         err = hypot(*x - x0, *y - y0);
60:         printf("%4d, % .15e, % .15e, % .15e\n", n, *x, *y, err);
61:     } while (err >= eps);
62:     return 0;
63: }
64: //  ------------    Main    ------------
65: main()
66: {
67:     double x=1, y=0.5, eps = 1e-10;
68:     char s[128];
69:
70:     fprintf(stderr, "# x0 y0 = "); gets(s); sscanf(s, "%lf%lf", &x, &y);
71:     Newton2(&x, &y, eps);
72:     printf("\n# (x y) = (% .15e % .15e)\n", x, y);
74:     return 0;
75: }

    "# x0 y0 = 1 0.5" と入力したときの実行結果とグラフは

    課題4-2

    上記2つのプログラムを以下の要求に沿ってできるところまで改良せよ。 満たせなかったものがあるときは、それをレポートに明示すること。 以下の(2)の2変数のニュートン法の課題をある程度できた人は、(1)については省略してよい。

    (1) 2.2 のニュートン法のプログラム Newton.c は、 2.1 の二分法や 2.3 の連立非線形方程式のプログラムのように 方程式の左辺の関数が独立していない。 そこで、以下の要求を満たすプログラムを作成せよ。
    (1-1) 方程式の左辺関数 f (x) (= x 2 - a ) および その導関数 df (x)/dx を独立させる。 ただし、入力変数 a はグローバル変数としてよい。
    (1-2) ニュートン法のアルゴリズムの部分 (20行目〜28行目)を Newton() と名づけ、main関数から独立させる。

    double Newton(double x, double eps)
{
    …………………………………
    ニュートン法のアルゴリズム
    …………………………………
    return x;
}
    (1-3) 方程式の左辺の関数 f (x) と その導関数の名前をmain関数側で自由に指定できるように、 Newton() の引数に関数ポインタを含めて定義する。
    double Newton(double x, double (*f)(double), double (*df)(double), double eps)
{
    …………………………………
    ニュートン法のアルゴリズム
    …………………………………
    return x;
}
    なお、実行例として a = 3 の場合について報告せよ。

    (2) 2.3 の2変数のニュートン法のプログラム Newton2.c は、 同じような関数が多くてソースが冗長である。 そこで、以下の要求を満たすプログラムを作成せよ。
    (2-1) 関数f1(), f2()を配列 f[1], f[2]、 関数J11()〜J22()を2次元配列 J[1][1]〜J[2][2] に割り当てることにして、 方程式の左辺関数f (x)と そのヤコビ行列 J (x)を おのおの一つのサブルーチン(func(), Jcalc())にまとめること。
    ただし、関数の戻り値で渡していた値は引数で受け渡すことにし、 手続きの戻り値は終了コード(整数)とする。 

    int func(double x, double y, double f[])
{
    ……………………
    f[1] = …………;
    f[2] = …………;
    return 0;
}
    また、Cでの配列の定義は0から始まることに注意し、一つ多めに宣言しておく。 配列の数(ND)をあらかじめ以下のように定数マクロで宣言しておくとよい。
    #define ND 3
    2次元配列を引数でうまく渡せなかった人は、 次善の策として1次元配列にして渡すか、グローバル変数で引き渡すこと。
    (2-2) 変数x, yを配列x[1], x[2]、 変数ux, uyをu[1], u[2]に割り当ててサブルーチン等の引数の数を減らし、 n 変数に移行しやすくする。
    (2-3) 方程式の左辺の関数 f (x) と その導関数の名前をmain関数側で自由に指定できるように、 Newton2() の引数に関数ポインタを含めて定義する。
    なお、実行例では
    　　x 2 + (y/2)2 = 1, 　　y 2 = x 2 + 1
    について、第1象限の解を求めたときの結果と途中経過のグラフを報告せよ。

    【目次】 | 【1.】 | 【2.のトップ】 | 【3.】 | 【付録1】 【付録2】