ニュートン法を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

＞正常には作動しているのですが
＞範囲がおかしいのか１回しか計算されない

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

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

a>0.0001

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

お礼

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

回答No.5

orcus0930

f(x) = (cos(x))^2 - 0.5の微分は、
f'(x) = -sin(2x)ですね。

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

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

お礼

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

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

回答No.4

R_Earl

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

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

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

となります。

お礼

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

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

回答No.3

ushioni

回答ではなく、Ｃの書き方のアドバイスですが、
#define f2(x) sin(2*x)
は
#define f2(x) sin(2*(x))
と書いたほうが良いです。

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

お礼

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

直してみます。

回答No.2

cametan_42

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がおっしゃられている書き方が基本だと理解してはいるのですが…。

回答No.1

koko_u_u

微分がまちがっている。

お礼

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

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

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