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

締切済

質問者：19830617
質問日時：2009/01/25 12:15
回答数：6件

方程式 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;

}

}

通報する

この質問への回答は締め切られました。

質問の本文を隠す

このQ&Aに関連する最新のQ&A

「C言語計算」に関するQ&A： C言語初心者で早速詰まってしまいました。 C言語の計算についてでちょっと良く分からないです… #in

さらに表示

「プログラム言語」に関するQ&A： Windows8のメモ帳に対応しているプログラム言語

「C言語ニュートン法」に関するQ&A：（C言語）ニュートン法のプログラムの為の数値微分

「ニュートン法初期値」に関するQ&A：ニュートン法について　初期値

回答 (6件)

最新から表示
回答順に表示

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
- 件

通報する

この回答へのお礼

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で確認できました。

- 0
- 件

通報する

この回答へのお礼

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

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

通報する

お礼日時：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

となります。

- 0
- 件

通報する

この回答へのお礼

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

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

通報する

お礼日時：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してみると分かります。

- 0
- 件

通報する

この回答へのお礼

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

直してみます。

通報する

お礼日時：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によるプログラムの場合、プログラムの記法自体にはさほど問題が見られません。つまり、こう言うケースでは「プログラムを記述した人間が」何かミスをした、って事ですね。
このソースも微分係数をキチンと書き換えれば、そのままキチンとした計算結果を返してくれました。つまり、全体としてのロジックは全然間違ってなく、単に「与えた微分係数が間違ってた」ってだけですね。

- 0
- 件

通報する

この回答へのお礼

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

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

通報する

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

No.1

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

微分がまちがっている。

- 0
- 件

通報する

この回答へのお礼

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

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

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

通報する

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

このQ&Aに関連する人気のQ&A

「ニュートン法初期値」に関するQ&A：仕事が続かない。働きたくない。

さらに表示

「プログラム言語」に関するQ&A： C言語のプログラムが実行できません。

「C言語ニュートン法」に関するQ&A：ニュートン法を使って解を求めるC言語プログラム

「C言語計算」に関するQ&A：Ｃ言語を独学でマスターした人。。

お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう！

質問する（無料）

このQ&Aを見た人はこんなQ&Aも見ています

専門家回答数ランキング

専門家

※過去一週間分の回答数ランキングです。

質問する（無料）

他の専門家の回答をチェック

このQ&Aを見た人が検索しているワード

このQ&Aと関連する良く見られている質問

Qニュートン法を使って解を求めるC言語プログラム

C言語を使って　y=x^2-4x　のｙの解をニュートン法を使って求める
プログラムを作る課題を出されたんですが、ニュートン法が良く分かっていないので、いろいろ調べたり、人に聞いたりしたところ
#include<stdio.h>
#include<math.h>
void main()
{
int counter=0;
double an,g,f,sh=0.0001;
printf("初期値を入力して下さい==>");
scanf("%ld",&an);
do{
g=(an*an)/(2*an-4);
f=2*an-4;
counter++;
}while(fabs(f)>sh);
printf("反復回数 %d 回　y=%lf \n",counter,g);
}
でプログラムがこんな感...続きを読む

Aベストアンサー

ニュートン法を理解するには、せめて微分の初歩の初歩は理解して下さい。
とりあえず、必要なところだけ。

y=x^2-4x
の解をもとめるには、次のようにします。
計算途中の近似解をa(n)とすると、それを元に得られる、より精度の良い近似解a(n+1)は、
a(n+1)=(a(n)*a(n))/(2*a(n)-4)
として求められます。（ニュートン法の理論から）
これを繰り返していって、
|a(n+1)-a(n)| < 0.0001
となったときに、十分精度の良い解が得られたと判断し、計算を終了します。
（計算終了の閾値0.0001は提示されたプログラムから取りました。）

プログラムの間違いは、下記の２点。
誤：scanf("%ld",&an);
正：scanf("%lf",&an);

誤：f=2*an-4;
正：f=g-an; an=g;

上記を修正し、初期値が2より大きい場合は4.000000が、初期値が2未満のときは
0.000000が求められることを確かめて下さい。

他の回答も見る

Q３次方程式の求解プログラム（ニュートン法）について！

初心者です！
ニュートン法での３次方程式の求解プログラムを作成しようとしています。
ですがよく分かりません(-_-;)
どなたか教えていただけないでしょうか。。。

Aベストアンサー

プログラム書きました。
これで、どうでしょう？

/* 3次方程式の求解 */
#include<stdio.h>
#include<math.h>
#include<float.h>

#define PI 3.14159265358979323846

#define REAL 0 /* 実根 */
#define IMAGE 1 /* 虚数根 */

int main(void){

double a, b, c, d; /* 3次方程式の係数 */
double p, q;
double D;
double x1, x2, x3; /* 根の実部 */
double y1, y2, y3;
double Im; /* 根の虚部 */
double u, v;
double d_u, d_v;
double phi;

int kind; /* 根の種類: REAL, IMAGE */

// ここから
printf("3次方程式(ax^3+bx^2+cx+d=0)の各係数：\n");
printf("a>>>");
scanf("%lf", &a);
printf("b>>>");
scanf("%lf", &b);
printf("c>>>");
scanf("%lf", &c);
printf("d>>>");
scanf("%lf", &d);

// ここまで消して、a, b, c, d値を代入する式を書いてください。

printf("3次方程式(%gx^3%+gx^2%+gx%+g=0)の解は\n", a, b, c, d );
if( a == 0.0 ){
printf("a(%g)がゼロなので，データエラーとみなし，\n", a );
printf("解を求めずに終了します．\n");
return 1; /* 終了コードを1にして，プログラムを終了 */
}

p = ( 3*a*c - b*b ) / ( 9*a*a );
q = ( 2*b*b*b - 9*a*b*c + 27*a*a*d ) / ( 54*a*a*a );

D = q*q + p*p*p; /* Dには誤差が含まれてしまうことに注意 */

if( fabs(D) <= DBL_EPSILON ){ /* Dを数値的にゼロかどうかチェックする */
/* 3つの重根(2重根・3重根)を含む実根 */
y1 = 2*pow(-q, 1.0/3.0);
y2 = y3 = -y1/2;
kind = REAL;
}else if( D < 0 ){ /* 3つの異なる実根 */
phi = acos( -q / pow( -p*p*p, 1.0/3.0 ) );
y1 = 2*sqrt(-p)*cos(phi/3);
y2 = -2*sqrt(-p)*cos(phi/3+PI/3);
y3 = -2*sqrt(-p)*cos(phi/3-PI/3);
kind = REAL;
}else{ /* 1つの実根，2つの共役複素数の解 */
d_u = -q+sqrt(D);
d_v = -q-sqrt(D);

if( d_u >= 0 ) u = pow( d_u, 1.0/3.0 );
else u = -pow( -d_u, 1.0/3.0 );
if( d_v >= 0 ) v = pow( d_v, 1.0/3.0 );
else v = -pow( -d_v, 1.0/3.0 );
y1 = u+v;
y2 = y3 = -y1/2;
Im = sqrt(3)/2*(u-v);
kind = IMAGE;
}

x1 = y1 - b/(3*a);
x2 = y2 - b/(3*a);
x3 = y3 - b/(3*a);

if( kind == REAL ){
printf("x1=%g, x2=%g, x3=%g\n", x1, x2, x3 );
}else{
printf("x1=%g, x2=(%g + %g i), x3=(%g - %g i)\n",
x1, x2, Im, x3, Im );
}

return 0;
}

他の回答も見る

Q三角関数の記述の仕方

タイトルそのまんまなんですが、三角関数はＣ言語ではどのように記述すればいいでしょうか？
角度にラジアン表記でπ（パイ）を使いたいんですが、その表記方法もわかりません。
僕の持っている本に載ってなかったので質問させていただきました。
よろしくお願いします。

Aベストアンサー

Ｃ言語で三角関数を使うためには、math.h をインクルードする必要があります。使い方は例えば、こんな感じです。

#define M_PI 3.14159265358979 /* 円周率 */

double x, y, theta;

theta = M_PI / 4.0;
x = cos(theta); /* sin,cos,tanの引数は弧度法の角度です。*/
y = sin(theta);

πは上記の例のように自分で定義して使ってください。

他の回答も見る

Q方程式を2分法を用いて解くプログラム

学校で出されたＣプログラムの課題で、１問だけどうしても出来ない問題があるんです。
「方程式 f(x) = x2 - 2 = 0 を 2 分法を用いて解くプログラムを作成せよ。ここで、方程式 f(x) = x2 - 2 は関数として定義せよ。上位の方から 4 桁目まで正しい値が出たらループを止めるようにする。」
というものなのですが、この「２分法」というやり方もよく分かりません。
プログラムの作成方法と併せて教えて頂けると幸いです。

Aベストアンサー

こんな感じでしょうか？

#include <stdio.h>

/* 関数 f(x) */
double f(double x) {
　return x*x-2.0;
}

/* 二分法初期値 x1<x2 と誤差限界 eps を入力 */
double bisec(double x1, double x2, double eps) {
　double x;
　while (x2 - x1 >= eps) {
　　x = (x1+x2)/2.0; /* 中点計算 */
　　if (f(x1)*f(x) > 0.0) { /* 同符号か判定 */
　　　x1 = x;
　　} else {
　　　x2 = x;
　　}
　}
　return (x1+x2)/2.0;
}

int main(void) {
　double eps=0.00001;
　printf("%lf %lf\n",bisec(-2,0,eps), bisec(0,2,eps));

　return 0;
}

他の回答も見る

QC言語のプログラムが実行できません。

C言語のプログラムが実行できません。

コンパイルは出来るんですが、実行すると、「Segmentation fault」と表示されてしまいます。

これは何のエラーなんでしょうか？
基本的な事かもしれませんが、分かる方宜しくお願い致します。

Aベストアンサー

僕も何度も出したなぁ。

ひとくちにSegmentation faultといっても、それこそさまざまな要因があるので、
これだけで原因を突き止めるのは非常に難しいです。

コンパイルはあくまで文法としてみているだけであり、
処理の流れ、メモリ確保など、プログラムそのものを見ているわけではありません。
このエラーが出るのは文法などよりもっと上位の原因なのです。
たとえばですが。
長さ10の配列があったとして、11番目以降を参照したりすると、
そういうのが出たような気がします。

ですから、変数があれば、その内容をprintf文で逐一出していき、
変な値が入っていないとか、少しずつ直していくしかないと思います。

他の回答も見る

QC言語の二分法のプログラムについて

二分法によりルート２の近似値を求めるプログラム、ってどうやって作ったらいいんですか？
ちなみに初期値は２で、収束条件は10^-5です。
収束までの回数も求めなきゃいけません(ニュートン法と比較するため。ちなみにニュートン法はできました。)

似たような質問を見つけたのですが、どれも、難しいプログラムばかりで解読ができません。
関数とかif else文とかwhile文とかfor文とか、そういう簡単なのしか習ってないので、それで作れる範囲で教えてくださる方、
いらっしゃいましたら、よろしくお願いします。

Aベストアンサー

これでどうでしょう。
if文判定、及びwhile文判定を修正しました。

#include <stdio.h>
#include <math.h>
#define f(x)(x*x-2.0)
int main(void)
{
int i=0;
double m,x1=2.0,x2=1.0,eps=1.0e-5;
do{
++i;
m=(x1+x2)/2;
if (f(m)>0)
{
x1 = m;
}
else
{
x2 = m;
}
}
while (fabs(x1-x2)>eps);
printf("%f,%d\n",m,i);
return 0;
}

他の回答も見る

Qプログラムでの数字につく”f”の意味

こんばんは。
プログラムで数字に"f"が付いている（4.0fなど)ものを見かけるのですが
この"f"の意味はなんなのでしょうか？

よろしくお願い致します。

Aベストアンサー

float型を意味します
一応float型は単精度/double型は倍精度なので、
コンパイラ（場合によっては同じコンパイラでも対象CPU毎に）によって
割り当てられるBit長(精度)が異なる場合があります。

他の回答も見る

QVisial C++おけるπの使い方

自宅でCプログラミングの練習をするためVisial C++ 2008を使って
プログラムをしています。y<sin(πx)となった時の
割合などを計算するプログラムで　　　　　　　　　　　　　　　　「M_PIが定義されていない識別子です」
とでてきます。所持している本を参考にしてもM_PI=πとして用いる
と書いてあり、math.hもインクルードしてるので原因が分からなくて
困ってます。Visial c++ではπの使い方には何か別の使い方がある
のでしょうか？よろしくお願いします。
*↓が実際に作ったプログラムです。
#include ...続きを読む

Aベストアンサー

★アドバイス
・math.hをインクルードする前に『_USE_MATH_DEFINES』定数を define します。

#define _USE_MATH_DEFINES
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

int main( void )
{
　：
　return 0;
}
必ず include する前に定義して下さい。

他の回答も見る

QC言語で-23乗を取り扱うには

C言語でボルツマン定数

1.38*10^(-23)

を取り扱いたいと思っています。

k = 1.38/pow(10, 23);
あるいは
k = 1.38/100000000000000000000000;

としても、コンパイルはできますが実行したときに
『1.#INF00』
という表示(マイナス無限大？）が出てしまいます。

このようにとても小さな定数をC言語で扱う場合、
どのようにすればよいのでしょうか？

よろしくお願いします。

Aベストアンサー

#include <stdio.h>

int main(void)
{
double k = 1.38e-23;

printf("ボルツマン定数K=%e\n", k);
return 0;
}

他の回答も見る

Q（C言語）ニュートン法のプログラムの為の数値微分

Visual StudioでC言語のプログラムを書く課題が出されたのですが、ニュートン法のプログラムを書くために
関数f(x)=3+ 2x + {2x - cos^2(x) - sin(x)} / {x^2 + cos(x^2) + sin^3(x^3)}
という関数を微分しなければいけません。担当の先生からは、関数が複雑で普通に微分するのは難しいから
f'(x) ≒ {f(x+Δx) - f(x)} / Δx を、Δx=0.0001として数値微分しなさいと指定されています。
どのように計算すればいいのでしょうか？
Δxにそのまま代入するだけではとても計算できそうにないです。

Aベストアンサー

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

//#define MAX_ITERATIONS 50 　/* ニュートン法の最大反復回数 */
#define MAXITERATIONS 50 //上だとエラーになるのでこう書き換えた

double function(double x);
double dFunction(double x);
double newtonSolver(double initValue, double tolerance);

void main(){
double x;
x = newtonSolver(4.0, 0.001); /* 1つ目の引数が初期値、2つ目は許容誤差 */
printf("function(%f) = %f\n", x, function(x)); // ０になるかのチェックだよ
}

/* ニュートン法による解の計算 */
/* initValue x の初期値 */
/* tolerane 許容誤差 */
double newtonSolver(double initValue, double tolerance){
double x = initValue; /* 初期値の代入 */
int counter; // C言語だとこうしないといけない

printf("starting Newton solver...\n");
printf("step\tx\t\tcorrection\t\n");

/* for (int counter = 0; counter < MAXITERATIONS; ++counter){ これはC++やJavaだ */
for (counter = 0; counter < MAXITERATIONS; ++counter){
/* 求める解への修正量を求める */
double correction = -function(x) / dFunction(x);
printf("%3d\t%5f\t%f\n", counter, x, correction);
x += correction;

/* 修正量(correction)の絶対値が許容誤差(tolerance)以下だったら
解となる x を返して関数を終了する
ちなみにfabs()は絶対値を求める関数 */
if (fabs(correction) < tolerance){
printf("ans = %f\n\n", x);
return x;
}
}
/* 最大反復回数を超えた場合は、収束していないので
異常値をあらわすNot a Number (NAN)を返す */
printf("Newton solver did not converge!\n\n");
return NAN;
}

/*対象の関数*/
double function(double x){
return 3 + 2*x + (2*x - pow(cos(x),2) - sin(x)) / (pow(x,2) + cos(pow(x,2)) + pow(pow(sin(x), 3),3));
}

//関数の部分はNo.3さんからのコピー（ポリポリ）
//間違っていたら、No.3さんのせいだからね～（笑）

/* 対象の関数を微分した関数 */
double dFunction(double x){
double dx;
dx = 0.0001;
return (function(x+dx)-function(x))/dx;
}

私が動かしたプログラムはこれ↑よ。
ダメならば、上のソースプログラムをコピペして、コンパイルしろ！！

計算結果は

starting Newton solver...
step x correction
0 4.000000 -6.135696
1 -2.135696 2.833077
2 0.697381 -1.283641
3 -0.586260 -0.494514
4 -1.080774 0.095784
5 -0.984990 0.014808
6 -0.970182 0.000271
ans = -0.969911

function(-0.969911) = -0.000000

です。