sinx+x=1の計算をしたい

Question

c言語でsinx+x=1の計算をするという課題が出たのですが、上手くいきません。以下がコードになっています。わかる方いらっしゃいましたらお答えください。

#include <stdio.h>
#include <math.h>
#include <float.h>
#pragma warning(disable: 4819)

int i;
double x = 0;
double ya = 0, yb = 0;

double main(void) {
	
	for (i=0; i < 1000000000; i++) {
		x += 0.000001;
		ya = sin(x);
		yb = 1-x;
		printf("%f,           %f\n", ya, yb);
		if (ya == yb) {
			printf("解は%fです", x);
			printf("その時のyは%fです", ya);
			break;
		}
	}
	return 0;
}

Comments

[mather]

「うまくいきません」という部分を具体的に「期待したことは何で、実際に起っていることは何で、試したことは何か」書いてください。
おそらく解が表示されないのだと思いますが…。

[mather]

なお、学校の課題に関することはまず担当教官や同じ授業を受けている人に質問をするようにしましょう。

Answer 1

sin(x) + x = 1 ⇔ xin(x) + x - 1 = 0 をニュートン法とかで解けばいいと思います。

理論的なしくみはニュートン法で調べれば解説記事がいくらでもでてくるので、そちらを参照ください。

ニュートン法とは何か？？ニュートン法で解く方程式の近似解 - Qiita

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

#define EPS 0.0000001
#define ITERATIONS 100

double f(double x)
{
    return sin(x) + x - 1;
}

double df(double x)
{
    return cos(x) + 1;
}

int main(void)
{
    int k;
    double x_new;
    double x = 0;  // 初期点
    printf("x0: %f", x);

    for (k = 0; k < ITERATIONS; ++k) {
        x_new = x - f(x) / df(x);  // 更新
        printf("%d: x=%f, x_new=%f\n", k, x, x_new);

        if (fabs(x - x_new) < fabs(x) * EPS)
            break;
        x = x_new;
    }

    if (k == ITERATIONS)
        printf("no solution\n");
    else
        printf("solution: x=%f\n", x);
}

x0: 0.0000000: x=0.000000, x_new=0.500000
1: x=0.500000, x_new=0.510958
2: x=0.510958, x_new=0.510973
3: x=0.510973, x_new=0.510973
solution: x=0.510973

追記

元のコードの ya == yb が以下の点で問題

1. 0.000001 ずつ +x 方向に移動して f(x) の値を確認するというやり方なので、完全に0になる点が見つかるとは限らない。
2. 浮動小数点数同士の演算では、丸め誤差が発生するので、等号を使った比較は意図した通りには動かない。十分近いかどうかで判断する必要がある。

浮動小数点数の誤差を考慮した比較【double/float型の正しい比較方法】 | MaryCore

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

#define EPS 0.00001
#define ITERATIONS 1000000000

int main()
{
    double x = 0;
    for (int i = 0; i < ITERATIONS; ++i) {
        double y = sin(x) + x - 1;

        if (fabs(y) < EPS) {
            printf("solution: f(%f) = %f\n", x, y);
            // solution: f(0.510969) = -0.000008
            break;
        }

        x += 0.000001;  // 移動
    }
    return 0;
}

Answer 2

.sin(x) + x = 1のような非線形方程式の数値解法には

1. 二分法
2. Newton法
3. 割線法
4. はさみうち法

などがあります。今回ははさみうち法で求めたいと思います。

Qiita 非線形方程式の数値解法

はさみうち法の特徴として

1. 微分係数を求める必要がない
2. Newton法の初期値を求めるのに適している

があります。

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

double f(double);
double calc(double, double);

int main(void){
	double alpha;
	alpha = calc(0,M_PI / 6.0);
	printf("%f \n", alpha);
	return 0;
}

double f(double x){
	return sin(x) + x - 1;
}

double calc(double a, double b){
	double c;
	do{
		c = ( a * f(b) - b * f(a) ) / ( f(b) - f(a) );
		if(f(c) == 0){
			break;
		}
		else if( f(a) * f(c) < 0){
			b = c;
		}
		else if( f(a) * f(c) > 0){
			a = c;
		}
	}while(fabs(f(c)) > 1e-10);
	
	return c;
}

結果　x = 0.510973

はさみうち法の特徴で「Newton法の初期値を求めるのに適している」があるので求めた解を利用して
ニュートン法で求めてみました。

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

double df(double);
double f(double);
double calc(double, double);
double newton(double);

int main(void){
	double alpha,alpha2;
	alpha = calc(0,M_PI / 6.0);
	alpha2 = newton(0);
	printf("はさみうち法 : %.20f \n", alpha);
	printf("ニュートン法 : %.20f \n", alpha2);
	return 0;
}

double f(double x){
	return sin(x) + x - 1;
}
double df(double x){
	return cos(x) + 1;
}

double calc(double a, double b){
	double c;
	do{
		c = ( a * f(b) - b * f(a) ) / ( f(b) - f(a) );
		if(f(c) == 0){
			break;
		}
		else if( f(a) * f(c) < 0){
			b = c;
		}
		else if( f(a) * f(c) > 0){
			a = c;
		}
	}while(fabs(f(c)) > 1e-10);
	
	return c;
}

double newton(double x){
	double new_x;
	while(1){
		new_x = x - f(x) / df(x);
		if(fabs( f(x) ) / df(x) < 1e-10){
			break;
		}
		if(fabs( f(x) ) < 1e-10){
			break;
		}
		x = new_x;
	}
	return new_x;
}

結果

はさみうち法 : 0.51097342939223233670
ニュートン法 : 0.51097342938856915580

Comments

[episteme]

"ニュートン"だろが"はさみうち"だろが、(強い)単純増加/減少関数じゃないと求まらんことがあるので注意。

Answer 3

「他の探索法でやれ」というのも回答でしょうけど，
とりあえず質問者様が行っている方法に近い処理方法も欲しいのではないかな，ということで，
【xを0から0.000001刻みで調べて行って「一番良いx」を答えにする】というものを．
こんな感じでしょうか．

int main(void)
{
    double BestX = 0; //一番良かったxの値を記録
    double BestAbsDiff = DBL_MAX;  //一番良かった際の | sin(x) - (1-x) | の値を記録

    double x = 0;
    while( 1 )
    {
        double ya = sin( x );
        double yb = 1.0 - x;

        double Diff = ya - yb;
        double AbsDiff = fabs( Diff );
        if( AbsDiff < BestAbsDiff )
        {//これまでよりも良いなら結果を更新
            BestX = x;
            BestAbsDiff = AbsDiff;
        }

        if( Diff >= 0.0 )break;  //(※)終了条件

        x += 0.000001;
    }
    printf( "解は%fです\n", BestX );
    printf( "その時のyは%fです\n", sin(BestX) );
    return 0;
}

(※)終了条件　に関しては

• y = sin(x)
• y = 1-x

の２本のグラフを描いて眺めてみれば，
xを0から増加させていくなら，sin(x)が1-xを上回ったら時点でもうこれ以上探さなくてよいことがわかるので，このようにしています．
（問題で扱う式自体が特定のものに指定されているなら，そういう知見(?)を実装しても良いかな，と）

Comments

[fana]

おっと，最後は sin(BestX) だけでなく， 1-BestX も表示すべきですね．