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

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

No.6 ベストアンサー 回答者： funaho 回答日時： 2006/07/17 11:48

プログラム書きました。

これで、どうでしょう？

/* 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;
}

この回答への補足

プログラムまで作っていただいてどうもすみません！
プログラムが↓の部分辺りからよくわからないのですが
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 );

あと、初期値も表示したいのですが(^_^;)
まだ、始めたばかりの初心者なので理解悪くてすみません！

補足日時：2006/07/17 15:04

No.7 回答者： funaho 回答日時： 2006/07/17 18:11

>プログラムが↓の部分辺りからよくわからないのですが

>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 );
ここは、()の中の3*a*cとb*bを先に計算してから引いてます。下も同じような感じです。
これは、演算子が優先順位が高い方から計算されるからです。

演算子の優先順位関して
http://www9.plala.or.jp/sgwr-t/c/sec14.html#s14-5

>あと、初期値も表示したいのですが(^_^;)
プログラム内に記述済み。
これです。「printf("3次方程式(%gx^3%+gx^2%+gx%+g=0)の解は\n", a, b, c, d );」

この回答へのお礼

どうもありがとうございました！
参考にして頑張ってみたいと思います。

お礼日時：2006/07/18 05:12

No.5 回答者： funaho 回答日時： 2006/07/16 08:28

>ax^3+bx^2+cx+d=0の場合についても解が求められ、初期値は入力しなくていいようにしたいのですが^_^;

入力をやめて、初期値を設定するようにすれば、問題ないと思います。

この回答への補足

わざわざありがとうございます。
a.b.c.dを入力しニュートン法で解を求める。初期値は入力しないという決まりで作成しなくてはいけないのです(>_<)

補足日時：2006/07/16 21:43

No.4 回答者： swimeye 回答日時： 2006/07/16 08:24

ニュートン法を使わずに解の公式から解析的に求めたほうが早いのでは？

３次方程式に限定するとどうなるかかんがえてませんが、
ニュートン法は収束性を保障しないので、
任意の式に対して、かならず解を得られるように
初期値を自動で決めることはできないような気がします

No.3 回答者： funaho 回答日時： 2006/07/15 20:40

これですか。

参考URL：http://www.tonami-h.tym.ed.jp/kadai/H17suu2.pdf

この回答への補足

ax^3+bx^2+cx+d=0の場合についても解が求められ、初期値は入力しなくていいようにしたいのですが^_^;
無理ですかね・・・？
何度もすみません<(_ _)>

補足日時：2006/07/16 06:44

No.2 回答者： funaho 回答日時： 2006/07/15 14:12

ここらへんを参考にしてみては。

参考URL：http://www.ma.noda.tus.ac.jp/u/tg/comp/L11.pdf

この回答へのお礼

わざわざありがとうございます<m(__)m>
現在、
１．f(x) = ax＾3 + bx＾2 + cx + dのa.b.c.dの値を入力する
２．判別式Dによって極値の数を求める
３．極値の数により場合分けをして初期値の位置を決める（初期値もプログラムが決定）
４．Newton法により解を求める
という考えでいます。
参考に載っていた数値計算プログラムの例とは少し違うように思うのですが、説明が難しくてなかなか理解できません(>_<)

お礼日時：2006/07/15 16:56

No.1 回答者： episteme 回答日時： 2006/07/15 12:14

なにが/どうわからないのかがわからないので答えようがありません。

この回答への補足

説明不足でした...すみません(>_<)
以下が私の考えです
１．f(x) = ax＾3 + bx＾2 + cx + dのa.b.c.dの値を入力する
２．判別式Dによって極値の数を求める
３．極値の数により場合分けをして初期値の位置を決める（初期値もプログラムが決定）
４．Newton法により解を求める
この３・４辺りが分からずに困っています。

補足日時：2006/07/15 17:00