ONLY DO WHAT ONLY YOU CAN DO

こけたら立ちなはれ 立ったら歩きなはれ

C++で非線形方程式を解く (はさみうち法)

非線形方程式の解法(はさみうち法)を利用して2の平方根を求める

f:id:fornext1119:20140531161845p:plain
考え方は、2分法とほとんど同じ.
1. まず, 条件 a < b, f(a) < 0, f(b) > 0 を満たす点 a, b を考えると,
関数 f(x) の解は, 区間 (a,b) の中に存在する.
2. 次に, 点 (a,f(a)) と 点 (b,f(b)) を結ぶ直線と x軸の交点 c を考えると,
f(c) < 0 であれば, 解は区間 (c,b) の中に存在し,
同様に, f(c) > 0 であれば, 区間 (a,c) の中に存在する.
3. この作業を繰り返して, 区間を絞り込んで行くことで解を求める.

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

double f(double x)
{
    return x * x - 2;  
}

int main()
{
    double low = 1;
    double high = 2;
    double fa = f(low);  
    double fb = f(high);  

    double x;
    while (true)
    {
        x = (low * f(high) - high * f(low)) / (f(high) - f(low));
        printf("%12.10f %12.10f %12.10f\n", high, low, x);
            
        double fx = f(x);  
        if (fabs(fx) < 0.0000000001) break;

        if ((fx > 0) == (fa > 0))
        {
            low = x;  
            fa = fx;  
        }
        else                  
        {  
            high = x;  
            fb = fx;  
        }
    }

    printf("%12.10f %12.10f\n", x, sqrt(2));
    return 0;
}
2.0000000000 1.0000000000 1.3333333333
2.0000000000 1.3333333333 1.4000000000
2.0000000000 1.4000000000 1.4117647059
2.0000000000 1.4117647059 1.4137931034
2.0000000000 1.4137931034 1.4141414141
2.0000000000 1.4141414141 1.4142011834
2.0000000000 1.4142011834 1.4142114385
2.0000000000 1.4142114385 1.4142131980
2.0000000000 1.4142131980 1.4142134999
2.0000000000 1.4142134999 1.4142135516
2.0000000000 1.4142135516 1.4142135605
2.0000000000 1.4142135605 1.4142135621
2.0000000000 1.4142135621 1.4142135623
2.0000000000 1.4142135623 1.4142135624
1.4142135624 1.4142135624

非線形方程式の解法(はさみうち法)を利用して2の立方根を求める

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

double f(double x)
{
    return x * x * x - 2;  
}

int main()
{
    double low = 1;
    double high = 2;
    double fa = f(low);  
    double fb = f(high);  

    double x;
    while (true)
    {
        x = (low * f(high) - high * f(low)) / (f(high) - f(low));
        printf("%12.10f %12.10f %12.10f\n", high, low, x);
            
        double fx = f(x);  
        if (fabs(fx) < 0.0000000001) break;

        if ((fx > 0) == (fa > 0))
        {
            low = x;  
            fa = fx;  
        }
        else                  
        {  
            high = x;  
            fb = fx;  
        }
    }

    printf("%12.10f %12.10f\n", x, cbrt(2));
    return 0;
}
2.0000000000 1.0000000000 1.1428571429
2.0000000000 1.1428571429 1.2096774194
2.0000000000 1.2096774194 1.2388370021
2.0000000000 1.2388370021 1.2511598712
2.0000000000 1.2511598712 1.2562955295
2.0000000000 1.2562955295 1.2584233391
2.0000000000 1.2584233391 1.2593027847
2.0000000000 1.2593027847 1.2596659013
2.0000000000 1.2596659013 1.2598157668
2.0000000000 1.2598157668 1.2598776087
2.0000000000 1.2598776087 1.2599031258
2.0000000000 1.2599031258 1.2599136544
2.0000000000 1.2599136544 1.2599179985
2.0000000000 1.2599179985 1.2599197909
2.0000000000 1.2599197909 1.2599205304
2.0000000000 1.2599205304 1.2599208356
2.0000000000 1.2599208356 1.2599209615
2.0000000000 1.2599209615 1.2599210134
2.0000000000 1.2599210134 1.2599210348
2.0000000000 1.2599210348 1.2599210437
2.0000000000 1.2599210437 1.2599210473
2.0000000000 1.2599210473 1.2599210488
2.0000000000 1.2599210488 1.2599210495
2.0000000000 1.2599210495 1.2599210497
2.0000000000 1.2599210497 1.2599210498
2.0000000000 1.2599210498 1.2599210499
2.0000000000 1.2599210499 1.2599210499
1.2599210499 1.2599210499
参考文献