京大理学部 物理科学課題演習 C1 （その２）

#2次方程式の解き方

　2次方程式には解の公式があるので別にプログラムを作るほどのこともないのですが、解を数値的に求めたいときはやはりプログラムを作って計算しなければなりません。2次方程式

$$ax^2 + bx + c = 0\ \ \ \ \ (a \neq 0)$$

の解は

$$x = \frac{-b \pm \sqrt{b^2 -4ac}}{2a}$$

ですね。まずはこの公式をそのまま使ってプログラムを書いてみましょう。方程式の係数はキーボードから入力することにします。

#include <stdio.h>
#include <math.h>

int main(void)
{
    float a, b, c;
    float D, x1, x2;

    printf("係数 a, b, c の入力: ");
    scanf("%f%f%f", &a, &b, &c);

    D = b * b - 4.0f * a * c;
    x1 = (- b + sqrtf(D)) / (2.0f * a);
    x2 = (- b - sqrtf(D)) / (2.0f * a);

    printf("%15.8e %15.8e\n", x1, x2);

    return 0;
}

　ここで sqrt() を使用しているので、数学関数のライブラリをリンクする必要があります。したがってコンパイル時に -lm オプションをつけなければなりません。ファイル名を niji.c とした場合

$ gcc niji.c -lm -o niji

としましょう。このプログラムを実行すると、係数の入力を促すプロンプトが出るので、係数の値を空白で区切って入力します。

　さて、このプログラムがすぐに破綻することは明らかです。虚数解を持つ場合を扱えません。次に虚数解も含めたプログラムを作るわけですが、その前にこのプログラムでひとつおもしろい実験をしてみましょう。（以下の説明はコメント欄で指摘されたとおり、桁落ちの説明としては不適当でした。どうぞコメント欄でのやり取りをご参照ください。）次のような方程式を考えます。

$$(x - 1.23456 \times 10^8)(x - 1.23456 \times 10^{-8}) = 0$$

すなわち

$$x^2 - 123456000.0000000123456x + 1.5241383936 = 0$$

です。とりあえず

$$a=1.0, \ \ \ b=- 1.234560000000000123456e8, \ \ \ c=1.5241383936$$
として実行してみましょう。結果はどうなりましたか？小さい方の解が正しく求まっていませんね。理由を考えてみましょう。まず大きい方の解は

$$x_1 = \frac{-b + \sqrt{b^2 -4ac}}{2a}$$

から計算されています。また小さい方の解は

$$x_2 = \frac{-b - \sqrt{b^2 -4ac}}{2a}$$

から計算されています。ここで右辺の分子に着目してみましょう。今の場合 $b$ は負ですから $x_1$ の右辺の分子は（正の数＋正の数）になっています。一方で $x_2$ の方は（正の数−正の数）になっています。ところで実験に使った方程式の係数だと、$-b$ と $\sqrt{b^2-4ac}$ がほとんど同じ値になっていますね（$b^2$ に比べ $4ac$ は無視できるぐらい小さい）。このように2つの数の値がほとんど同じ場合、それらの間で引き算をすると、有効桁数が著しく失われてしまいます。例えば $9.87654321$ から $9.87654312$ を引くと、結果は $0.00000009$ となります。はじめ9桁の有効桁を持っていた数が、互いに引き算をした結果、有効桁がわずか1桁になってしまいました。これでは計算の精度が保てません。こういった理由で、小さい方の解が正確に計算できなかったわけです。このことを桁落ちといいます。

　では、桁落ちを避けるにはどうしたらよいでしょうか。引き算を含む方の解（今の場合は $x_2$ の方ですが、当然これは $b$ の符号によります）を別な方法で求めればよいわけです。ここでは解と係数の関係を使って $x_1$ から

$$x_2 = c\ /a\ /x_1$$

として求めることにしましょう。$b$ の符号が正の場合は先に $x_2$ を求め、$x_1$ を上と同様に $x_2$ から求めます。このことを考慮したプログラムを以下に示します。

#include <stdio.h>
#include <math.h>

int main(void)
{
    float a, b, c;
    float D, x1, x2;

    printf("係数 a, b, c の入力: ");
    scanf("%f%f%f", &a, &b, &c);

    D = b * b - 4.0f * a * c;
    if (b >= 0.0) {
        x2 = (- b - sqrtf(D)) / (2.0f * a);
        x1 = c / a / x2;
    }
    else /*if (b < 0.0)*/ {
        x1 = (- b + sqrtf(D)) / (2.0f * a);
        x2 = c / a / x1;
    }

    printf("%15.8e %15.8e\n", x1, x2);

    return 0;
}

　プログラムを入力して実行してみてください。今度は両方の解が正確に求まるはずです。

　さて、いよいよ虚数解を含む完全なプログラムを作る段階にきました。虚数解を扱うためには、判別式

$$D = b^2 - 4ac$$

の値によって処理の制御をしなければなりません。以下にプログラム例を示します。各自入力して、係数を変えて試してみてください。

#include <stdio.h>
#include <math.h>

int main(void)
{
    double a, b, c;
    double D;
    double x1, x2;
    double xrprt, xiprt;

    printf("係数 a, b, c の入力: ");
    scanf("%lf%lf%lf", &a, &b, &c);

    D = b * b - 4.0 * a * c;

    if (D > 0.0) {
        if (b >= 0.0) {
            x2 = (- b - sqrt(D)) / (2.0 * a);
            x1 = c / a / x2;
        }
        else /*if (b < 0.0)*/ {
            x1 = (- b + sqrt(D)) / (2.0 * a);
            x2 = c / a / x1;
        }
        printf("相異なる実数解\n");
        printf("%15.8e %15.8e\n", x1, x2);
    }
    else if (D == 0.0) {
        printf("重解\n");
        printf("%15.8e\n", - b / (2.0 * a));
    }
    else if (D < 0.0) {
        xrprt = -b / (2.0 * a);
        xiprt = sqrt(-D) / (2.0 * a);
        printf("虚数解\n");
        printf("%15.8e %15.8e i\n", xrprt, xiprt);
        printf("%15.8e %15.8e i\n", xrprt, -xiprt);
    }

    return 0;
}

　上記のプログラムには、実は決定的な欠陥があります。それは、入力値として $a=0$ を読み込んだ場合に正しく計算ができないことです。そこで、 $a=0$ を読み込んだときには、エラーメッセージを出力してプログラムが停止するように修正してあげる必要があります。

　さて前回、ひとまとまりの手続きは関数にするとよい、と力説しました。そこで、ここでは2次方程式を解く関数を作ることにしましょう。関数の名前は qdrteq とします。まずはじめに、この関数の機能を定めておきます。この関数は、2次方程式の係数（$a, b, c$）を受け取り、判別式の値（$D$）と方程式の解（$x_1, x_2$）を返すことにします。さらに関数の終了状態を表す変数 icond も返すことにします。すなわち、正常に解を求められた場合は icond = 0 を返し、係数に誤りがあり（$a=0$）解が求められなかった場合は icond = 1 を返すようにします。このように、エラーの処理も関数に行わせるわけです。係数の入力と結果の表示はメイン関数で行うことにします。当然エラーがあった場合の表示もメイン関数で行います。解の受け渡しの方法ですが、ここでは2次元の配列を2つ（ x_1, x_2 ）用意し、それぞれの配列の第一成分に解の実部を、第二成分に虚部を代入することにします。もちろん解が実数解であった場合は、虚部には $0$ がはいることになります。これでプログラムの仕様がだいたい決まりました。

##課題

　上の解説を参考にして、2次方程式の解を求める関数を完成させなさい。また、2次方程式の係数を入力する部分と、結果を表示する部分も関数にしなさい。

　以上が実現されると、メイン関数には変数の宣言と関数を呼び出す文しか残らないことになります。係数を入力し、解を求め、結果を表示するという手続きの流れのみがメイン関数に明確に記述されるわけです（手続きの具体的な内容は各関数に書かれている）。これで最高に読みやすいプログラムが完成します。