第5回非線形方程式の解法

5.1 まずは、解くべき方程式を把握する

非線形方程式の解を解析的に求めることは非常に困難であるが、計算機を用いて数値的に求める方法がいくつか知られている。
非線形方程式の解を求める場合、多くは反復計算を行う。従って収束判定が必要となる。解くべき方程式を、f(x) = 0 とする。関数 f(x) が、区間 a ≦ x ≦ b でどのような値を取るかを出力せよ。

アプローチ

関数 f(x) = 2x³ -4x + 1 を対象の関数とする。以下のように、関数 f(x) をC言語の関数として定義し、利用すること。

double func(double x) {
	return( 2*x*x*x - 4*x + 1 );
}

ここでは、対象とする区間を、a = -2, b = 2として考える。このとき、値の推移をなるべく細かく把握できるよう、この区間を n=100 分割して値を表示することにする。

#include<stdio.h>

double func(double x) { ... (上に書いてあるので、省略)

int main(void) {
  const double a = -2.0;
  const double b =  2.0;
  const int n = 100;

  for(int i=0; i<=n; i++) { ...
(以降の処理を考えてください)

以下のような出力結果が得られるよう、上のプログラムの続きを作成せよ。

x = -2.000, f(x) =   -7.000
x = -1.960, f(x) =   -6.219
x = -1.920, f(x) =   -5.476
...
x =  1.960, f(x) =    8.219
x =  2.000, f(x) =    9.000

結果を確認

表示結果だけでも、求める解は大まかに把握することができるが、ついでにExcelでグラフを書いてみましょう。
まずは、以上の結果をファイルに書き出すことにする。もちろん fprintf を利用しても良いが、今回はリダイレクトと呼ばれる機能を利用して、画面に出力される結果をファイルに書き出してみる。(良い手抜き)

z:\compmech\ex04> func_test.exe > data.txt

これらの値を利用すると、Excelなどを利用して下図のようにグラフの形状を確認することもできる。 f(x) = 0 を満たす解は、区間[-2 -1], [0 1], [1 2]に存在することがわかる。

図解くべき方程式の把握

5.2 関数を引数とする関数

これまでは関数の引数として、数値やポインタを利用してきたが、多くのプログラミング言語では、「関数を引数とした関数」を定義することができる(高階関数という)。
以下のようなプログラムで、高階関数を利用した例である。このプログラムでは、main文の show_func(func, -2, 2, 100); という命令で、関数 func を引数とした関数 show_func を利用していることに注意せよ。
# 実は、関数を引数とする関数も、結局はポインタを渡していることになる。

#include<stdio.h>

typedef double fp(double);	//関数ポインタ

double func(double x) {
	return( 2*x*x*x - 4*x + 1 );
}

double func2(double x) { 
	return( 4*x*x*x - 3*x*x + 2*x + 1 );
}

void show_func(fp* func, const double a, const double b, const int n) {
  for(int i=0; i<=n; i++) {  ...
  (5.1と同じ処理を記述)
  }
}

int main(void) {
  show_func(func, -2, 2, 100);
  show_func(func2, -2, 2, 100);
  return 0;
}

5.3 二分法による非線形方程式の解法

f(x) = 0 を満たす解を数値的に求める方法の一つとして、二分法が知られている。
二分法は以下のように行う。

図二分法

二分法のアルゴリズム

解の左右にある2点a,bをlow,highの初期値とする。a,bの区間はなるべく狭く設定しておくこと。
lowとhighの中点xをx=(low+high)/2で求める。
f(x)>0 なら解はxより左にあるからhigh=xとし、上限を半分に狭める。
f(x)<0 なら解はxより右にあるからlow=xとし、下限を半分に狭める。
f(x)が0または|high-low|< eps となったときのxを求める解とする。
そうで無ければ 2. へ戻って計算を続ける。
※このときの eps を収束判定値という。例えば eps = 10^(-6) などと十分に小さな値としておく。

これをプログラムで書くと、次のようになる。

プログラム例 (PDFファイル)

5.4 二分法のプログラムの修正

二分法を用いて解を求めるためには、f(a)>0かつf(b)<0 の場合も含めて考慮する必要がある。上のサンプルプログラムでは、f(a)<0かつf(b)>0 であるとしていた。
f(a), f(x), f(b) の符号関係に着目し、f(a)<0かつf(b)>0 の場合と、f(a)>0かつf(b)<0 の場合の、どちらの場合にも利用できるプログラムを作成すること。

5.5 ニュートン法による非線形方程式の解法

ニュートン法(ニュートン・ラプソン法とも呼ばれる)では以下の手順で解を求める。

図ニュートン法

ニュートン法のアルゴリズム

解の近くの適当な値を初期値とする．
のにおける接線を引き，軸と交わったところをとし，以下同様の手順でと求めていく．
- ここで、からを求める方法は、点、傾きの直線の方程式を考えることにより
  
  である。
になったときのの値を求める解とする．そうでなければ2.以降を繰り返す．

課題 (提出〆切：金曜日22時)

以下の関数において、-2.0 ≦ x ≦ 2.0 の範囲で方程式 f(x) = 0 を満たす全ての解を画面に出力するプログラムを作成せよ。

以下の指示に従うこと。
- ニュートン法により解を求めること。
- 全ての解に対して、初期推定値と反復回数も表示すること。
- 初期推定値は、コマンドライン引数から与える形式にすること。
- プログラム名は 5_1.cpp とする。
- 求める関数の導関数はあらかじめユーザ側で準備しておく必要がある。この課題では以下のような関数を定義して利用すること。
```
// 関数f(x)の微分式
double dfunc(double x)
{
  return /* 求める関数の導関数 */ ;
}
```
- 以下は実行例である。値は適当なので注意すること。
```
Z:\compmech\ch05>5_1.exe 0.12 0.45 0.78

----- Answers ‘f(x) = 0.0’ -----
x[1] = 1.234		(initial value: 0.12  -->  iteration: 3 times)
x[2] = 5.678		(initial value: 0.45  -->  iteration: 6 times)
x[3] = 9.012		(initial value: 0.78  -->  iteration: 9 times)
```
図に示すようなスライダクランク機構がある。スライダの位置 x とクランクの回転角θの関係は、次式で示される。

r = 1.0, l = 2.0 のとき、0 ≦ θ ≦ π の範囲で、x = 1.5 となる θ を求めるプログラムを作成せよ。以下の指示に従うこと。
- 2分法とニュートン法の 2 通りの方法で解を求めること。
- 2 分法における解の範囲の上限と下限、およびニュートン法における初期推定値はコマンドライン引数から与える形式にすること。解の下限、解の上限、初期推定値の順に引数を指定すること。
- 解に対して、解の範囲の上限と下限、初期推定値、反復回数を表示すること。
- プログラム名は 5_2.cpp とする。
- 以下は実行例である。値は適当なので注意すること。
```
Z:\compmech\ch05>5_2.exe 1.2 3.4 1.2

----- Answers ‘x = 1.5’ -----
< Bisection Method >
theta = 1.234 rad	(low: 1.2, high: 3.4  -->  iteration: 56 times)
< Newton Method >
theta = 1.234 rad	(initial value: 1.2  -->  iteration: 3 times)
```

第5回 非線形方程式の解法