ニュートン法をC言語でプログラム

方程式 cos^2x-0.5=0 (0<x<π) の解をニュートン法で求める
という問題をC言語のプログラムを作り計算したいのですが分かりません。
自分で考えてみたプログラムは以下の通りです。

#include <stdio.h>
#include <math.h>
#define f1(x) cos(x)*cos(x)-0.5
#define f2(x) sin(2*x)

/* ニュートン法による方程式の解 */

main()
{
double x0,x1,a,b,c,d,g,n;

a=1;
x0=0.7;
n=0;

while(a>0.0001){
b=x0;
d=f1(b);
g=f2(b);
x1=x0-d/g;
c=x1;
a=f1(c);
n=n+1;

printf(" n= %f x1=%f x0=%f\n",n,x1,x0);
printf(" a= %f → 解 x= %f \n", a,x1);

x0=x1;

}

}

自分としてはこれが精一杯で、何故間違ってるのか、何をどうすればいいのか、さっぱり分かりません。どういったところが間違ってるのか可能性だけでも示して頂ければ幸いです。

参考として、ニュートン法によるプログラム例として書かれていたものを上げさせて頂きます。
例： e^x-3=0 の解をニュートン法により計算する。

#include <stdio.h>
#include <math.h>
#define f1(x) exp(x)-3
#define f2(x) exp(x)

/* ニュートン法による方程式の解 */

main()
{
double x0,x1,e,a,b,c,d,g,n;

a=1;
x0=3;
n=0;

while(a>0.0001){
b=x0;
d=f1(b);
g=f2(b);
x1=x0-d/g;
c=x1;
a=f1(c);
n=n+1;
printf(" n= %f x1=%f x0=%f\n",n,x1,x0);
printf(" a= %f →　解　　x= %f \n", a,x1);

x0=x1;

}

}

No.6 回答者： cametan_42 回答日時： 2009/01/25 20:00

＞正常には作動しているのですが

＞範囲がおかしいのか１回しか計算されない

そうですか?
初期値x0に解の近傍の0.7与えているんで1回しか計算されなくてもしょーがない、とか思うんですけどね。
一般に、ニュートン法は収束までのステップ数は短いんで、当然だ、と思いますよ。

それより、問題は収束判定の

a>0.0001

の方だと思います。これがaが負になった場合、対処してない、って思うんですが。
例えば、計算でf1(c)＝-2000何とやら、となった場合、そこで計算が止まってしまいます。これは要求仕様じゃないでしょ?
多分、やりたかったのはaの絶対値>0.0001だと踏んでるんですが、それがコードに反映されていません。問題の性質から言うと、初期値はπの付近(例えば3.0とか)からスタートしても解を返さなきゃならない筈、なんですが、元々のコードを微分の部分だけ手直ししたら、おかしな解が出てくるんじゃないのかな、って思います。
その辺、ちょっと見てみてください。

この回答へのお礼

0.7だから１回しかされないのと絶対値にするの分かりました。ありがとうございました。

お礼日時：2009/02/01 02:03

No.5 回答者： orcus0930 回答日時： 2009/01/25 16:34

f(x) = (cos(x))^2 - 0.5の微分は、

f'(x) = -sin(2x)ですね。

あと、引数の代入で、double型に代入するときは
a=1;
ではなく、
a=1.0;
とするほうが無難でしょう。

あなたのプログラムをちょっと修正するだけで、正常作動することは、私のPCで確認できました。

この回答へのお礼

回答ありがとうございました。

正常には作動しているのですが
範囲がおかしいのか１回しか計算されないので
プログラムで計算させる意味がない事になってしまいました。

お礼日時：2009/01/25 18:01

No.4 回答者： R_Earl 回答日時： 2009/01/25 16:23

> (sin^2x)'

> =2(sinx)(sinx)'
> =2sinx*cosx
> =sin(2X)
>
> となると思ったのですが…
> よろしければ微分の何が違ってるかご指摘下さい。お願いします。

微分される関数はcos^2x - 0.5ですよね？
なぜsin^2xの微分をしているのでしょうか？

(cos^2x - 0.5)'
= (cosx)(-sinx)
= -(sinx)(cosx)
= -sin2x

となります。

この回答へのお礼

何度もありがとうございます。
焦って違う問題の計算メモを書いてしまいました。
すみません。
マイナスがなかったのも自分の字が汚かっただけのようです。

ご丁寧にありがとうございました。

お礼日時：2009/01/25 17:56

No.3 回答者： ushioni 回答日時： 2009/01/25 13:12

回答ではなく、Ｃの書き方のアドバイスですが、

#define f2(x) sin(2*x)
は
#define f2(x) sin(2*(x))
と書いたほうが良いです。

理由は、
f2(PI/2)とf2(0+PI/2)をprintfしてみると分かります。

この回答へのお礼

ご指摘ありがとうございます。

直してみます。

お礼日時：2009/01/25 14:48

No.2 回答者： cametan_42 回答日時： 2009/01/25 12:49

Cのマナーの部分と、数学の部分と、2つありますね。

Cの部分:

1:書き出しはmain()じゃなくってint main(void)にした方が良い。

これはツマラない部分ですし、計算結果は同じですが、「Cのマナー」としては、引数を取らないmain関数では引数はvoid、とした方が良い、です。
また、main関数の返り値自体も厳密に言うと、「計算結果とは別にして」OSに対して整数が返されます。従って、

2:正常にプログラムが終了した事を伝える為、返り値を返すように必ず最後にreturn 0;を付ける。

これもマナーです。計算結果を別にOSに「正常終了したよ～～。」と教える為、return 0;も必ず記述するようにしましょう。

数学の部分:

#1さんも仰っている通り、(cos(x))^2-0.5の微分はsin(2*x)になりません。計算ミス、ですね。
プログラムは「言われた通りに計算してる」だけですんで、与えた微分係数が間違ってても感知せずに、「言われた通りに計算している」だけです。
特に、コンパイル時にエラーが出ないで正常実行出来るCによるプログラムの場合、プログラムの記法自体にはさほど問題が見られません。つまり、こう言うケースでは「プログラムを記述した人間が」何かミスをした、って事ですね。
このソースも微分係数をキチンと書き換えれば、そのままキチンとした計算結果を返してくれました。つまり、全体としてのロジックは全然間違ってなく、単に「与えた微分係数が間違ってた」ってだけですね。

この回答へのお礼

回答ありがとうございます。

Cのマナーについてですが、今このように授業で習っているのでその先生に対して違う書き方でプログラムを書いて出すことは出来ないのでどうかそこは見逃して下さい。一応、cametan_42がおっしゃられている書き方が基本だと理解してはいるのですが…。

お礼日時：2009/01/25 14:40

No.1 回答者： koko_u_u 回答日時： 2009/01/25 12:34

微分がまちがっている。

この回答へのお礼

さっそく回答ありがとうございます

(sin^2x)'
=2(sinx)(sinx)'
=2sinx*cosx
=sin(2X)

となると思ったのですが…
よろしければ微分の何が違ってるかご指摘下さい。お願いします。

お礼日時：2009/01/25 14:34