連立方程式を解きたい
■ このスレッドは過去ログ倉庫に格納されています
0001大学院生
NGNG連立方程式をN×Nの行列にしてから、
ガウスの消去法かガウス・ジョルダン法を使おうと思っていたのですが、
除算を使うため、誤差が出てしまいます。
なにか良い方法はありませんか？
開発環境はVCです。
0002名無しさん＠お腹いっぱい。
NGNGそうでないのなら、体をクラスを作ればいいんじゃない？。浮動小数点を
よくわからんが、、、
体をクラスにしたらどうよ？
例えば、有理数を係数とする連立方程式なら、
0003名無しさん＠お腹いっぱい。
NGNG0004２
NGNG要は、数学でいうところの「体」をクラスにしたらどうか、ということ。
0005名無しさん＠お腹いっぱい。
NGNGそれか他のゼミ生に相談するとか。
0006大学院生
NGNG厳密って言うのは誤差がないっていうことです。
体をクラスにっていうのは、すいません、ちょっとわかりません。
3の方の除算を分数のままというのは最初考えたのですが、
複雑になりすぎてしまうような気がするのですが、
やるしかないのでしょうか？
数値計算の本は一応一通り読んでみました。
しかし、どの本も厳密ということはあまり意識していないようでした。
0007名無しさん＠お腹いっぱい。
NGNG問題は以下に許容範囲ないに誤差を納めるか。
0008大学院生
NGNGそうなんですけど、最終的な答えを分数で出せば、
最小限に押さえられるかなと。
でも、やっぱり3さんのいうように、全て分数で行うしかないのかな。
0009名無しさん＠お腹いっぱい。
NGNGMathematicaは確か任意の精度で実数演算できたような気がします。
それに行列演算もサポートしてますし。
ただ速度は劇遅です。
0010名無しさん＠お腹いっぱい。
NGNG0011２
NGNGどうやら、有理数を係数とする方程式を解くことが目的のようだな。
順番として、
(1)double型で誤差があってもいいから連立方程式を解くコードをかけ。
(2)有理数体クラス（分数クラス）を定義しろ。
連立方程式を解くのに必要となる operator（=＠｀ + ＠｀-＠｀ *＠｀ / など）
を作れ。
(3)(1)で作ったコードの、doubleを有理数体クラスに置き換えろ。
以上。
0012Visual名無しさん
NGNGそれはともかく、もし連立方程式の解を知ること自体が目的なら>>9の
言うように、MathematicaなりMapleなり使うべきだと思うよ。
0013名無しさん＠お腹いっぱい。
NGNG0014名無しさん＠お腹いっぱい。
NGNG0015名無しさん＠お腹いっぱい。
NGNG本質的には何も変わらないような...
0016名無しさん＠お腹いっぱい。
NGNGサイズがちょっと大きくなると破綻しそうだが
疎な行列だったらある程度単純化できそうな気がする。
0017名無しさん＠お腹いっぱい。
NGNGこれで破綻する大きさなら掃き出し等である程度の
小行列に対角化したりしておくという方法もあると思う。
0018名無しさん＠お腹いっぱい。
NGNGいずれにせよ１の扱う行列の癖というかクラスというか
使えそうな性質があったら書いた方がいいと思う。
0019名無しさん＠お腹いっぱい。
NGNGそのへんも書いたほうがいいんじゃないかな。
あと、どのくらいの時間でもとめるか？とかね。
厳密解がどうしても欲しければ、>>3 のとおり、
多倍精度整数の組で表される分数でごりごり計算するんだね。遅そうだけど。
そういう数学計算用のパッケージも世の中にごろごろあるし探してきて使えば？
個人的な興味では、精度保証付の浮動小数点演算をするのが面白そうだとおもう。
0020名無しさん＠お腹いっぱい。
NGNG小数を含む連立方程式の場合、浮動小数点数で扱ったらその時点で誤差が出るような？
0.1とかって２進数では正確に表現できないはず。
文字列で受け取って小数の場合は小数点以下の桁数を調べて分数表現に変換しないとマズイ。（1／10みたいな）
というより整数10の場合も 10／1 とか分数表現に統一して、分母を変数A、分子を整数B
みたいにそれぞれ変数を割り当てて計算すれば簡単に出きるのでは？
0021学歴っておいしい？
NGNG#include <stdlib.h>
#include <float.h>
#include <conio.h>
#define SUC 1
#define ERR 0
#define MAX 10
#define MIN 1
#define mabs( x ) ( (x > 0) ? x : -(x) )
void main( void );
char get_zisuu( char max＠｀ char min );
double *get_matrix( char zisuu );
double *get_vecter( char zisuu );
void set_matrix( char zisuu＠｀ double *mat );
void set_vecter( char zisuu＠｀ double *vec );
void put_vecter( char zisuu＠｀ double *vec );
void mswap( double *＠｀ double * );
char pratial( char zisuu＠｀ char gyou＠｀ double *matrix＠｀ double *vecter );
char gauss_j( char zisuu＠｀ double *matrix＠｀ double *vecter ); void main()
{
char zisuu;
double *matrix＠｀ *vecter;
while( 1 )
{
zisuu = get_zisuu( MAX＠｀ MIN );
if( ( matrix = get_matrix( zisuu ) ) == (double *)NULL ||
( vecter = get_vecter( zisuu ) ) == (double *)NULL )
{
puts( " メモリ−が取得できません。" );
}
else
{
set_matrix( zisuu＠｀ matrix );
set_vecter( zisuu＠｀ vecter );
if( gauss_j( zisuu＠｀ matrix＠｀ vecter ) == ERR )
{
puts( " 入力値エラ−です。" );
}
else
{
printf( " 解は>>>" );
put_vecter( zisuu＠｀ vecter );
}
}
free( (void *)matrix );
free( (void *)vecter );
}
}
char get_zisuu( char max＠｀ char min )
{
char zisuu;
while( 1 )
{
printf( " %d から %d までの数字で次数を入力して下さい、終了は 0 です。\n"＠｀ min＠｀ max );
if( scanf( "%d"＠｀ &zisuu ) != 1 )
{
scanf( "%*s" );
puts( " 数字を入力して下さい。" );
continue;
}
if( (min > zisuu) || (zisuu > max) )
{
if( !zisuu )
{
puts( "END" );
exit( 1 );/* END */
}
puts( " 決められた範囲内で入力して下さい。" );
continue;
}
return zisuu;
}
}
0022学歴っておいしい？
NGNG{
return (double *)malloc( sizeof( double ) * zisuu * zisuu );
}
double *get_vecter( char zisuu )
{
return (double *)malloc( sizeof( double ) * zisuu );
}
void set_matrix( char zisuu＠｀ double *mat )
{
char gyou＠｀ retu;
for( gyou = 0; gyou < zisuu; gyou++ )
{
for( retu = 0; retu < zisuu; retu++ )
{
while( 1 )
{
printf( " 行列 %d 行 %d 列の要素は？>>>"＠｀ gyou＠｀ retu );
if( scanf( "%lf"＠｀ mat + (gyou * zisuu) + retu ) == 1 )
{
break;
}
puts( " 入力ミスです。" );
scanf( "%*s" );
}
}
}
}
void set_vecter( char zisuu＠｀ double *vec )
{
char gyou;
for( gyou = 0; gyou < zisuu; gyou++ )
{
while( 1 )
{
printf( " ベクトルの %d 行の要素は？>>>"＠｀ gyou );
if( scanf( "%lf"＠｀ vec + gyou ) == 1 )
{
break;
}
puts( " 入力ミスです。" );
scanf( "%*s" );
}
}
}
void put_vecter( char zisuu＠｀ double *vec )
{
char gyou;
for( gyou = 0; gyou < zisuu; gyou++ )
{
printf( " %lf"＠｀ *(vec + gyou) );
}
puts( "" );
}
void mswap( double *a＠｀ double *b )
{
double damy;
damy = *a;
*a = *b;
*b = damy;
}
char partial( char zisuu＠｀ char gyou＠｀ double *matrix＠｀ double *vecter )
{
char roop＠｀ r_max;
double pivot;
roop = r_max = gyou;
pivot = *( matrix + (zisuu * roop) + gyou );
while( roop < zisuu )
{
if( mabs( pivot ) < mabs( *( matrix + (zisuu * roop) + gyou ) ) )
{
r_max = roop;
pivot = *( matrix + (zisuu * roop) + gyou );
}
roop++;
}
if( mabs( pivot ) < DBL_EPSILON )
{
return ERR;
}
if( gyou != r_max )
{
for( roop = 0; roop < zisuu; roop++ )
{
mswap( matrix + (zisuu * gyou) + roop＠｀ matrix +(zisuu * r_max) + roop );
}
mswap( vecter + gyou＠｀ vecter + r_max );
}
return SUC;
}
0023学歴っておいしい？
NGNG{
char gyou＠｀ roop＠｀ ro_2;
double pivot＠｀ damy;
for( gyou = 0; gyou < zisuu; gyou++ )
{
if( partial( zisuu＠｀ gyou＠｀ matrix＠｀ vecter ) == ERR )
{
return ERR;
}
pivot = *( matrix + (zisuu * gyou) + gyou );
for( roop = 0; roop < zisuu; roop++ )
{
*( matrix + (zisuu * gyou) + roop ) /= pivot;
}
*( vecter + gyou ) /= pivot;
for( roop = 0; roop < zisuu; roop++ )
{
if( roop != gyou )
{
damy = *( matrix + (zisuu * roop) + gyou );
for( ro_2 = 0; ro_2 < zisuu; ro_2++ )
{
*( matrix + (zisuu * roop) + ro_2 ) -= damy * *( matrix + (zisuu * gyou) + ro_2 );
}
*( vecter + roop ) -= damy * *( vecter + gyou );
}
}
}
return SUC;
}char gauss( char zisuu＠｀ double *matrix＠｀ double *vecter )
{
char gyou;
char damy;
char roop＠｀ ro_2;
for( gyou = 0; gyou < zisuu; gyou++ )
{
if( partial( zisuu＠｀ gyou＠｀ matrix＠｀ vecter ) == ERR )
{
return ERR;
}
for( roop = gyou + 1; roop < zisuu; roop++ )
{
damy = *( matrix + ( zisuu * roop ) + gyou ) / *( matrix + ( zisuu * gyou ) + gyou );
for( ro_2 = gyou; ro_2 < zisuu; ro_2++ )
{
*( matrix + ( zisuu * roop ) + ro_2 ) -= damy * *( matrix + ( zisuu * gyou ) + ro_2 );
}
*( vecter + roop ) -= damy * *( vecter + gyou );
}
}
}
0025かんそう
NGNG綴りは正しく書きましょう(roop＠｀damy＠｀vecter)
サイズの指定にcharを使うのはやめましょう。intかsize_tを使いましょう。
数学ライブラリの中に入出力のコードを書くのはやめましょう。
NULLにキャストは不要です。
mainの戻り値はvoidではありません。
0026名無しさん＠お腹いっぱい。
NGNG>サイズの指定にcharを使うのはやめましょう。intかsize_tを使いましょう。
>数学ライブラリの中に入出力のコードを書くのはやめましょう。
この２つは正しい。
あとは単なる煽り！
0027名無しさん＠お腹いっぱい。
NGNG行を交換してるから、順番を元に戻す処理を入れた方が親切かもしれない。
（意味取り違えてたらゴメン）
0028ため息
NGNGこれは煽りというより、万が一>>21-23 がどこぞの大学院生だとしたら
土下座して親に謝罪すべきというほうが正解だね。
loop＠｀ vector dummy程度を満足に綴れない大学生なんて、
はっきりいってゴミ。いったい何をやってきたんだい？
ところでそのダメコードの添削が希望なの？
0029名無しさん＠お腹いっぱい。
NGNGあまりにも、あんまりなんで、てっきりネタだと思ったよ。
0030名無しさん＠お腹いっぱい。
NGNG複素数クラスが倍精度実数２つを引数にとるように、64ビットの整数２つを
引数にとるクラスで、四則演算を定義する（複素数の有理数なら整数４ついるのか）。
注意としては演算や代入の度に約分していかないとすぐに桁が一杯になる。
でもそんな大きなＮは無理だと思うな（桁あふれするから）。まあ10×10くらいが
関の山では？
0031>28
NGNGこの程度のが多いよ。マジ逝って欲しい。
0032名無しさん＠お腹いっぱい。
NGNGある必要があるんですか? 数論ならともかく、どんな値にでも
多かれ少なかれ誤差があるので、連立方程式自体が無誤差で解けても
問題を無誤差で解いたことにはならないのでは?
既出の精度保証付き数値計算 (ISBN4-535-78253-X)とかのほうが
現実には使えるのではないかと思います。
>30 いちど有理数計算クラス作って、演算ごとに約分してましたが、
すんごい遅くて泣きました。約分できない場合も多く簡単にケタ溢れ
するんで、多倍長整数クラスの実装も必須になると思います。
新規情報少ないのでsage
0033名無しさん＠お腹いっぱい。
NGNG0034２１
NGNG勉強になるのでもっともんく書いて下さい
0035大学院生
NGNG研究内容はあんまり詳しくはかけないのですが、点列を色々な曲線で補間して、その滑らかさを比較する、
ということをしています。点列を曲線で補間するときに連立方程式を使います。
自分の実験用にしか使わないので、速度は遅くてもぜんぜん問題はありません。
やはり、有理数計算クラスを作るしかないみたいですね。
頑張って勉強してみます。
学歴っておいしい？さん、コードを書いていただいてありがとうございます。
参考にしてみます。
0036名無しさん＠お腹いっぱい。
NGNG>ということをしています。点列を曲線で補間するときに
>連立方程式を使います。
ふむふむ。数学的な議論で終始するなら、連立方程式の厳密評価でも
イイかもしんないですね。
でももしアルゴリズムの研究なら「そのアルゴリズムがどのくらい
簡単に使えるか」も、アルゴリズムの評価の一端になると思います。
「多くの場合、非常に特異に近い係数行列になりますが、
その連立方程式が厳密に評価できれば良い補間です」とかだと、
実用上は使い物になりませんよね。
一般的な(陰的ピボット選択付きのLU分解などの)連立方程式
解法による評価も併せて行うべきだと思います。
それから、新規の手法の評価に有理数計算クラスをつくるなら、
徹底的にテストしないと、手法に不備があるのか、計算クラスに
不備があるのかで、わけわかんなくなるので要注意。
がんばってください。
0037名無しさん＠お腹いっぱい。
NGNGhttp://www.interq.or.jp/jazz/iijima/equations/equations00.html
■ このスレッドは過去ログ倉庫に格納されています