C言語関係掲示板

過去ログ

No７２７　超越関数 (2^x) - tan(x) = 0　の解を求める

[戻る]　[ホームページ]
No.8918
関数の計算です
投稿者---ryi(2003/08/13 14:14:04)


超越関数 (2^x) - tan(x) = 0　の解を求めるプログラムがわかりません。わかる人がいたら教えてください。
No.8919
Re:関数の計算です
投稿者---たいちう(2003/08/13 15:14:15)


解を探すのは微分可能な部分だろうしNewton法が簡単でしょう。
もし知らなきゃ調べてください。
No.8925
Re:関数の計算です
投稿者---7of9(2003/08/14 02:02:17)


>解を探すのは微分可能な部分だろうしNewton法が簡単でしょう。
>もし知らなきゃ調べてください。

求まるのは、近似解です。

＃解を求めるプログラム？？？
＃多分近似解でいいんだろうな、^^;
No.8928
Re:関数の計算です
投稿者---かずま(2003/08/14 03:03:22)


> 超越関数 (2^x) - tan(x) = 0　の解を求めるプログラムがわかりません。

y = 2^x のグラフと、y = tan(x) のグラムは、簡単に書けますよね。
さて、y = tan(x) は、区間 (-π/2, π/2) で連続です。
f(x) = 2^x - tan(x) とすると、上記の 2つのグラフの交点の x座標が解です。
f(x) は、区間 (-π/2, π/2) で連続であり、かつ単調増加であることは
グラフを見ればすぐわかりますから、この区間での解を求めてみましょう。

Newton法を使うと収束が速いのですが、二分探索でも解を求めることが出来ます。
こうすれば、f(x) を微分する必要がありません。2^10 = 1024、2^20 = 1メガ、
2^30 = 1ギガですから、30回もループを回れば、10進で 9桁の精度は出るでしょう。
#include <stdio.h>
#include <math.h>

#define N    50
#define EPS  1e-9

double f(double x) { return pow(2, x) - tan(x); }

int main(void)
{
    int i;  double a = -1.57, b = 1.57, x, y;

    if (f(a) * f(b) > 0) return 1;
    for (i = 0; i < N; i++) {
        x = (a + b) / 2;
        y = f(x);
        if (fabs(y) < EPS) break;
        if (f(a) * y > 0) a = x; else b = x;
        printf("%3d: %12.9f\n", i, x);
    }
    return 0;
}

No.8930
Re:関数の計算です
投稿者---かずま(2003/08/14 09:07:06)


> f(x) = 2^x - tan(x) とすると、上記の 2つのグラフの交点の x座標が解です。
> f(x) は、区間 (-π/2, π/2) で連続であり、かつ単調増加であることは
> グラフを見ればすぐわかりますから、この区間での解を求めてみましょう。

ちょっと訂正。
f(x) は、2^x と tan(x) の差ですから、「単調増加」ではなく「単調減少」ですね。
No.8957
Re:関数の計算です
投稿者---ryi(2003/08/16 15:37:40)


if (f(a) * f(b) > 0) return 1;


if (fabs(y) < EPS) break;

上の二つの文の意味、必要性がよくわからないのですが
No.8958
Re:関数の計算です
投稿者---ryi(2003/08/16 15:55:47)


あと if (f(a) * y > 0)
この条件がよくわかりません　なぜf(a)*yなのか
No.8960
Re:関数の計算です
投稿者---かずま(2003/08/16 22:51:08)


> if (f(a) * f(b) > 0) return 1;

これは、f(a) と f(b) の符号が異なることを確認するためのもので、
この場合は不要でしょう。

区間 (-π/2, π/2) に f(x) = 0 の解があるといっても、
f(-π/2) や f(π/2) の値は存在しませんから、
区間 [-1.57, 1.57] で解を探そうとしました。
ひょとすると、区間 (1.57, π/2) に解が存在するかも知れず、
その場合、f(-1.57) と f(1.57) が同符号になるので、それを
検出するためにこれを書いてしまいました。


> if (fabs(y) < EPS) break;

これは、間違いですね。

区間 [a, b] をどんどん二分して、f(x) = 0 になる x に近づけて
行くのですが、ループを回っている途中で、たまたま y = 0 になって
しまえば、解が求まってしまったのですから、ループを回る必要は
なくなります。 そこで、if (y == 0) break; と書こうとしたのですが、
浮動小数点数の比較で == になることは稀なので、1e-9 程度の差なら
0 に等しいとしていいだろうと思って、if (fabs(y) < EPS) break;
と書いてしまいました。

でも、これでは、x が真の値に近づく前にループを抜けてしまうことも
ありますね。if (y == 0) break; にしてください。


> あと if (f(a) * y > 0)
> この条件がよくわかりません　なぜf(a)*yなのか

これは、二分探索で必須の条件です。
f(a) と f(b) は異符号です。
x = (a+b)/2 の f(x) がどちらの符号と同じかで、次の区間を
[a, x] にするのか、[x, b] にするのかが変ってきます。


それから、ループの回数ですが、if (y == 0) break; で抜けることは
普通ないでしょうから、もっとちゃんと指定しなければなりません。

30回まわると、2^30 = 1ギガですから、区間 [-1.57, 1.57] の幅 3.14
の 10億分の 1 の精度になります。

あるいは、ループの回数を適当に指定して、ループの中で、
if (b - a < EPS) break; を実行して、実際に区間の幅が十分小さく
なったことを検出してループを抜けてもよいでしょう。
No.8962
Re:関数の計算です
投稿者---かずま(2003/08/17 12:13:07)


> Newton法を使うと収束が速いのですが、二分探索でも解を求めることが出来ます。

Newton法だと、f(x) の微分 f'(x) = log(2) * 2^x - (sec(x))^2 を
使いますが、初期値 x0 の与え方で、収束が遅くなったり、求まる
解の区間が変ったりしますね。
#include <stdio.h>
#include <math.h>

double f(double x)  { return pow(2, x) - tan(x); }
double f1(double x) { return log(2)*pow(2, x) - 1/(cos(x)*cos(x)); }

#define EPS  1e-9

int main(void)
{
    double x0, x;

    for (x0 = 1; ; x0 = x) {
        x = x0 - f(x0)/f1(x0);
        if (fabs(x - x0) < EPS) break;
    }
    printf("%12.9f\n", x);
    return 0;
}