C言語の3次方程式

[east end girl]

A.
#include <stdio.h>

main()
{
float xa,xb,xc;
float ya,yb,yc;
float initial_delta,dx;


printf(" Xa Xc Xb 誤差 Ya Yc Yb\n");

initial_delta=1.0;
xa = 0.0;
xb = xa + initial_delta;
dx = initial_delta / 2.0;
xc = xa + dx;

ya = xa*xa*xa - 0.5;
yb = xb*xb*xb - 0.5;
yc = xc*xc*xc - 0.5;
printf("%10.7f %10.7f %10.7f %10.7f %10.7f %10.7f %10.7f\n",xa,xc,xb,dx,ya,yc,yb);
/* ----- ここから ----- */
if ( ( ya * yc ) < 0 ) {
xb = xc;
} else {
xa = xc;
}
dx /= 2.0;
xc = xa + dx;

ya = xa*xa*xa - 0.5;
yb = xb*xb*xb - 0.5;
yc = xc*xc*xc - 0.5;
printf("%10.7f %10.7f %10.7f %10.7f %10.7f %10.7f %10.7f\n",xa,xc,xb,dx,ya,yc,yb);
/* ----- ここまでをコピーする ----- */
}

B.Aで作ったプログラムを，while あるいは for 文による繰り返し処理を使って，同じ文がなるべくプログラム中に重複して出てこないように修正し，解を誤差 0.0001 以内で求めよ．その際，同じ意味を持つ定数がなるべく重複して出てこないようにする等，わかり易く効率的なプログラムにする修正も含めるものとする(スペースを適当に入れるだけでもわかり易さは格段に向上する)．出力の形式は自由で良いが，見やすいものにすること．さらに，プログラムにはある程度のコメントを入れること．元から入っているコメントは削除して構わない．


C.Bで作ったプログラムを参考にし，「0.83X^3+4.55x^2+5.19x+2.45=0」の実数解を，誤差 0.00002 以内で求めるプログラムを作れ(保存するディレクトリは ~/win_home 以下の任意のディレクトリとし，ファイル名は各自で適当に決めること)．出力の形式は自由で良いが，見やすい出力とすること．プログラムにはある程度のコメントを入れること．




という大学の課題が出たのですが、さっぱりです。C言語の授業５回目のレポートで１０時間粘ってみたのですが、何にも解決出来なかったため、ここで質問させて頂きます。どなたか助言をして頂けると助かります。

[みけCAT]

コードはcodeタグで囲んでいただき、適切にインデントをしていただくと、見やすくてありがたいです。

まず、このコードはどのような方程式の計算をするものですか？

[みけCAT]

ヒント
・for 文による繰り返し処理（問題文より）
・配列
・double型（誤差が少なくなる）

「0.83X^3+4.55x^2+5.19x+2.45=0」のX(大文字)とx(小文字)は別の変数とみなしますか？

[east end girl]

返信ありがとうございます！
この方程式は「x^3-0.5=0』及び「0.83x^3+4.55x^2+5.19x+2.45=0」の実数解を以下のステップを踏むことにおり、数値的に求めよと書いてありました。

for文による繰り返し処理で、結構苦戦しています。どこにforを置くべきかすらもわからなくなってきました。
また、変数は大文字小文字で別の変数とみなしません、ごめんなさい。

[みけCAT]

Aのプログラムは与えられたものですか？それとも自分で書いたものですか？

コード:

ya = xa*xa*xa - 0.5;
yb = xb*xb*xb - 0.5;
yc = xc*xc*xc - 0.5;

のように、同じような文が並んでいる場所がありますね。
x,yを配列にし、これらをfor文で処理するとよさそうです。

[east end girl]

Aのプログラムは与えられたものです！

for ( init_exp ; cond_exp ; update_exp ) {
statement1;
statement2;
......

for ( init_exp ; cond_exp ; update_exp )
statement;

↑これらを真似して行ったのですが、ﾁｮｯﾄ難しいですね・・・。

[みけCAT]

三次方程式の解 - 高精度計算サイトによると、x^3-0.5=0の解は

コード:

x = -0.39685026299205 +0.6873648184993i
x = -0.39685026299205 -0.6873648184993i
x = 0.7937005259841

の3個です。
また、Aのプログラムの出力は、Ideoneによると

コード:

Xa Xc Xb 誤差 Ya Yc Yb
 0.0000000  0.5000000  1.0000000  0.5000000 -0.5000000 -0.3750000  0.5000000
 0.5000000  0.7500000  1.0000000  0.2500000 -0.3750000 -0.0781250  0.5000000

です。
出力の意味を教えていただけますか？

[みけCAT]

for ( init_exp ; cond_exp ; update_exp ) {
statement1;
statement2;
......

for ( init_exp ; cond_exp ; update_exp )
statement;

↑これらを真似して行ったのですが、ﾁｮｯﾄ難しいですね・・・。

何を行ったか、教えていただけますか？

[east end girl]

printf(" Xa Xc Xb 誤差 Ya Yc Yb\n");の前にforを置いただけです・・・。

[NorthDipper]

３次方程式の解の公式を調べたところ、とても複雑な式で求めていたためその公式は使ってないと思われ、コードを見る限り２分探索のような方法で求めているものと思います。
この場合は方程式上の座標（xc,yc）においてyc＝０となるときのx座標であるxcが解になります（条件を満たすxcは複数ある場合が多い）。
このxcを求めるために、yaとybの正負が逆であり　xa<xc かつ xc<xb という条件を満たす方程式上の二点(xa,ya),(xb,yb)を定めると、yc=0 となるようなxcが必ずxaとxbの間に必ずあることになります。
あとは２分探索の要領でxaとxbの距離を縮めていくといずれ方程式の解となるxcが求まります。

[みけCAT]

NorthDipperさんが解説している方法には、欠点があるかもしれません。
二分探索は単調増加/減少するものに使うアルゴリズムですが、
xaとxbの値によっては図のような状況になり、この条件を満たさないことが考えられます。

三次関数の例

sanzi_test.png (9.31 KiB) 閲覧数: 8928 回

(雑な図ですみません)

そこで、対象の方程式を一度微分し、導関数を求めます。これは二次関数になるので、値が0になるxが簡単に求められます。
その情報を用いて二分探索の範囲を決めれば、この方法で解けるかもしれません。

導関数=0となるxをx1,x2(x1<x2)とおくと
・「x1の左側でx=x1のときのもとの関数の値と符号が逆の場所(x1-1,x1-1-2,x1-1-2-4, ...と調べていくとよさそうかも)」～x1
・x1～x2
・x2～「x2の右側でx=x2のときのもとの関数の値と符号が逆の場所(x2+1,x2+1+2,x2+1+2+4, ...と調べていくとよさそうかも)」
の三つの範囲で二分探索をすれば、解が出そうだと思います。
あとは、コーナーケースの処理を適宜してください。

[usao]

>実数解を以下のステップを踏むことにおり、数値的に求めよと書いてありました。
とのことですから，あれこれ想像で回答するよりも
まずはその指定されたステップとやらを確認された方が早い気が…

課題文（の一部だと思われる文面）を見た感じだと，課題としてやらせたいことは非常に初歩的なことで，
数値計算アルゴリズムの部分まで求めていない可能性もありそうですし．
（たとえばプログラムAのように，「ちゃんと答えに収束するような初期状態を与える」
といった類の前提条件があるのかもしれませんし．）

[NorthDipper]

みけCATさんの言う通り関数の極値の座標を求める手順を飛ばしていました。
関数の形によっては解が３つ出てこない場合もあると思うので、導関数を使って極値の座標を求めるときに関数の増減も求めて配列に記憶しておくとさまざまな形の３次方程式に対応しやすいと思います。