ニュートン法

ケプラー方程式x-e*sin(x)-c=0の解をステップ数とともに出力するプログラムで、e,cはそれぞれ0.5と1です。
xに値を入力して計算させるのですが、どうしてもできません。
下のプログラムリストでおかしいところはどこでしょうか？
// ニュートン法 x-e*sin(x)-c=0
#include <stdio.h>
#include <math.h>
#define e 0.5
#define c 1.0
#define K 10000

double fun(double x);
double bibun(int i,double x);
float m=1.0,n=1.0;
int i=1;

main(){
float x1,x2;
float z;
printf("初期値x0を入力して下さい\n");
scanf("%f",&x1);
for(i=0;i<=K;i++){
x2 = x1 - fun(x1)/bibun(i,x1);
x1 = x2;
z = fun(x1);
z = fabs(z);
if(fabs(z)<=0.00001){
break;
}
if(i==K){
printf("収束しません\n");
exit(0);
}
　　 }
printf("解 = %f\n",x1);
printf("ステップ数 = %d",i);
return 0;
}

// 関数f(x)
double fun(double x1){
double r;

r = x1 - e * sin(x1) - c;
return r;
}

// 微分
double bibun(int i,double x1){
float p;

if(i%2==1){
p = pow((-1.0),m)*e*sin(x1);
m++;
} else {
p = pow((-1.0),n)*e*cos(x1);
n++;
}
return p;
}

No.3 回答者： koko_u 回答日時： 2007/02/10 00:30

>koko_uさんの考えてるアルゴリズムはどのようなものですか？

アルゴリズムはニュートン法なんでしょ。

どうもわかっておられないようなので、ひとまず初期値を 0 とか置いて、 x_0 = 0, x_1 = ◯, x_2 = △, ... と計算して下さい。

その後、計算が面倒なのでコンピュータにやらせるとしたらどうすればいいか考えましょう。

この回答へのお礼

わかりました。
新ｘについて、新f(x)を微分するということをやっていました。
ニュートン法自体をうまく把握できていなかったようです。
ありがとうございました。

お礼日時：2007/02/10 22:10

No.2 回答者： noname#26650 回答日時： 2007/02/09 09:54

サンプルです。

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

#define e (0.5)
#define c (1.0)
#define K (10000)
#define EPS (0.00001)

double f(double x);
double g(double x);

int main(void)
{
　　double x1, x2;
　　int i;
　　
　　printf("初期値x1を入力して下さい\n");
　　scanf("%lf", &x1);
　　for (i = 0; i < K; i++) {
　　　　x2 = x1 - f(x1) / g(x1);
　　　　if (fabs(x2 - x1) <= EPS)
　　　　　　break;
　　　　x1 = x2;
　　}
　　if (i < K) {
　　　　printf("解 = %f\n", x1);
　　　　printf("ステップ数 = %d\n", i);
　　}
　　else {
　　　　printf("収束しません\n");
　　　　exit(1);
　　}
　　return 0;
}

double f(double x)
{
　　return x - e * sin(x) - c;
}

double g(double x)
{
　　return 1 - e * cos(x);
}


（注）インデントのため、全角空白を使っています。

No.1 回答者： koko_u 回答日時： 2007/02/09 00:44

とりあえず、double bibun(int i, double x1) がまったくわからん

1. 戻り値が double なのに何故か float p;
2. m, n が唐突に出現。
3. (-1)^m を求めようとして pow( -1.0, m ) としてるの？普通に場合分けじゃね？

だいたい、ニュートン法って一階微分がわかれば x_(n+1) = x_n - f(x_n)/f'(x_n) でステップバイステップじゃないの？
もっとややこしいことしようとしてるの？

この回答への補足

x_(n+1) = x_n - f(x_n)/f'(x_n) でステップバイステップというのがよく分かりません。x_(n+1)を新たなx_nと考えて計算していくんですよね？そうするとsinとcosが交互に出てくると思って、こんな感じにしてしまったのですが。
koko_uさんの考えてるアルゴリズムはどのようなものですか？

補足日時：2007/02/09 02:01