■掲示板に戻る■ ■過去ログ倉庫めにゅーに戻る■

連立方程式を解きたい

1 名前： 大学院生 投稿日： 2000/10/25(水) 19:39

研究で厳密な連立方程式の解を求めることになりました。
連立方程式をN×Nの行列にしてから、
ガウスの消去法かガウス・ジョルダン法を使おうと思っていたのですが、
除算を使うため、誤差が出てしまいます。
なにか良い方法はありませんか？
開発環境はVCです。

2 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/25(水) 19:59

厳密の意味がよくわからんが、実数（浮動小数点）を使うということか？

そうでないのなら、体をクラスを作ればいいんじゃない？。浮動小数点を
よくわからんが、、、
体をクラスにしたらどうよ？
例えば、有理数を係数とする連立方程式なら、

3 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/25(水) 20:04

除算を分数のままひきずり廻すとか。

4 名前： ２ 投稿日： 2000/10/25(水) 20:06

スマソ。作りかけの文を書いてしまった。
要は、数学でいうところの「体」をクラスにしたらどうか、ということ。

5 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/25(水) 20:22

素直に数値計算の本読んだ方がよいと思うが。
それか他のゼミ生に相談するとか。

6 名前： 大学院生 投稿日： 2000/10/25(水) 20:48

double型だと数値が丸められて誤差が出てしまうような気がするのです。
厳密って言うのは誤差がないっていうことです。
体をクラスにっていうのは、すいません、ちょっとわかりません。
3の方の除算を分数のままというのは最初考えたのですが、
複雑になりすぎてしまうような気がするのですが、
やるしかないのでしょうか？
数値計算の本は一応一通り読んでみました。
しかし、どの本も厳密ということはあまり意識していないようでした。

7 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/25(水) 20:57

コンピュータで計算するかぎり誤差ナシなんて無理だろ。
問題は以下に許容範囲ないに誤差を納めるか。

8 名前： 大学院生 投稿日： 2000/10/25(水) 21:01

>7
そうなんですけど、最終的な答えを分数で出せば、
最小限に押さえられるかなと。
でも、やっぱり3さんのいうように、全て分数で行うしかないのかな。

9 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/25(水) 21:04

VCじゃないとだめですか？
Mathematicaは確か任意の精度で実数演算できたような気がします。
それに行列演算もサポートしてますし。
ただ速度は劇遅です。

10 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/25(水) 21:12

Octaveはどのくらいの精度なんだっけ？

11 名前： ２ 投稿日： 2000/10/25(水) 21:17

「体をクラスに」ていうのは、３と同じこと。
どうやら、有理数を係数とする方程式を解くことが目的のようだな。

順番として、
(1)double型で誤差があってもいいから連立方程式を解くコードをかけ。
(2)有理数体クラス（分数クラス）を定義しろ。
　　連立方程式を解くのに必要となる operator（=, + ,-, *, / など）
　　を作れ。
(3)(1)で作ったコードの、doubleを有理数体クラスに置き換えろ。

以上。

12 名前： Visual名無しさん 投稿日： 2000/10/25(水) 21:51

「体」がわからん大学院生って困りもんだなぁ…

それはともかく、もし連立方程式の解を知ること自体が目的なら>>9の
言うように、MathematicaなりMapleなり使うべきだと思うよ。

13 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/25(水) 21:58

最大値軸選択や、スケーリングとかじゃあだめ？

14 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/25(水) 22:03

除算の誤差が気になるなら多倍長演算はどうよ?

15 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/25(水) 22:07

>14
本質的には何も変わらないような...

16 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/25(水) 22:24

余因子行列を定義通り計算するのは？

サイズがちょっと大きくなると破綻しそうだが
疎な行列だったらある程度単純化できそうな気がする。

17 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/25(水) 22:45

俺は数値計算の勉強したことないんで適当だが、
これで破綻する大きさなら掃き出し等である程度の
小行列に対角化したりしておくという方法もあると思う。

18 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/25(水) 23:04

「対角化」までは要らんな。ミス。

いずれにせよ１の扱う行列の癖というかクラスというか
使えそうな性質があったら書いた方がいいと思う。

19 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/25(水) 23:45

>>18 のとおり行列の大きさ(N)や性質(疎かどうか？)によるから、
そのへんも書いたほうがいいんじゃないかな。

あと、どのくらいの時間でもとめるか？とかね。
厳密解がどうしても欲しければ、>>3 のとおり、
多倍精度整数の組で表される分数でごりごり計算するんだね。遅そうだけど。
そういう数学計算用のパッケージも世の中にごろごろあるし探してきて使えば？

個人的な興味では、精度保証付の浮動小数点演算をするのが面白そうだとおもう。

20 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/26(木) 01:56

>>1
小数を含む連立方程式の場合、浮動小数点数で扱ったらその時点で誤差が出るような？
0.1とかって２進数では正確に表現できないはず。
文字列で受け取って小数の場合は小数点以下の桁数を調べて分数表現に変換しないとマズイ。（1／10みたいな）
というより整数10の場合も 10／1 とか分数表現に統一して、分母を変数A、分子を整数B
みたいにそれぞれ変数を割り当てて計算すれば簡単に出きるのでは？

21 名前： 学歴っておいしい？ 投稿日： 2000/10/26(木) 02:11

include <stdio.h>
#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;
}
}

22 名前： 学歴っておいしい？ 投稿日： 2000/10/26(木) 02:12

double *get_matrix( char zisuu )
{
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;

}

23 名前： 学歴っておいしい？ 投稿日： 2000/10/26(木) 02:12

char gauss_j( char zisuu, double *matrix, double *vecter )
{
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 );
}
}
}

24 名前： 学歴っておいしい？ 投稿日： 2000/10/26(木) 02:14

>>23　　　　　ガウス法とガウス所留弾法のプログラム
>>22,21　　　２３を操作するプログラム　　

25 名前： かんそう 投稿日： 2000/10/26(木) 02:21

識別子にローマ字を使うのはやめましょう（zisuu,gyou）
綴りは正しく書きましょう(roop,damy,vecter)
サイズの指定にcharを使うのはやめましょう。intかsize_tを使いましょう。
数学ライブラリの中に入出力のコードを書くのはやめましょう。
NULLにキャストは不要です。
mainの戻り値はvoidではありません。

26 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/26(木) 02:26

>>25
>サイズの指定にcharを使うのはやめましょう。intかsize_tを使いましょう。
>数学ライブラリの中に入出力のコードを書くのはやめましょう。
この２つは正しい。

あとは単なる煽り！

27 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/26(木) 04:39

> 22
行を交換してるから、順番を元に戻す処理を入れた方が親切かもしれない。
（意味取り違えてたらゴメン）

28 名前： ため息 投稿日： 2000/10/26(木) 06:42

＞綴りは正しく書きましょう(roop,damy,vecter)
これは煽りというより、万が一>>21-23 がどこぞの大学院生だとしたら
土下座して親に謝罪すべきというほうが正解だね。
loop, vector dummy程度を満足に綴れない大学生なんて、
はっきりいってゴミ。いったい何をやってきたんだい？
ところでそのダメコードの添削が希望なの？

29 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/26(木) 08:35

＜２４
あまりにも、あんまりなんで、てっきりネタだと思ったよ。

30 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/26(木) 09:50

まあ前に出てるようなことだけど、有理数計算クラスを作るのがいいのかな。
複素数クラスが倍精度実数２つを引数にとるように、64ビットの整数２つを
引数にとるクラスで、四則演算を定義する（複素数の有理数なら整数４ついるのか）。
注意としては演算や代入の度に約分していかないとすぐに桁が一杯になる。
でもそんな大きなＮは無理だと思うな（桁あふれするから）。まあ10×10くらいが
関の山では？

31 名前： >28 投稿日： 2000/10/26(木) 09:56

いや、数学科・物理学科の教授のコードだって
この程度のが多いよ。マジ逝って欲しい。

32 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/26(木) 11:41

>1 研究の内容がわからんのでなんとも言えませんが、本当に厳密で
ある必要があるんですか? 数論ならともかく、どんな値にでも
多かれ少なかれ誤差があるので、連立方程式自体が無誤差で解けても
問題を無誤差で解いたことにはならないのでは?
既出の精度保証付き数値計算 (ISBN4-535-78253-X)とかのほうが
現実には使えるのではないかと思います。

>30 いちど有理数計算クラス作って、演算ごとに約分してましたが、
すんごい遅くて泣きました。約分できない場合も多く簡単にケタ溢れ
するんで、多倍長整数クラスの実装も必須になると思います。

新規情報少ないのでsage

33 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/26(木) 18:10

２１〜２３以外で誰かこ−ど書いて下さい

34 名前： ２１ 投稿日： 2000/10/26(木) 18:14

大学１年生です
勉強になるのでもっともんく書いて下さい

35 名前： 大学院生 投稿日： 2000/10/26(木) 18:20

皆さん情報ありがとうございました。
研究内容はあんまり詳しくはかけないのですが、点列を色々な曲線で補間して、その滑らかさを比較する、
ということをしています。点列を曲線で補間するときに連立方程式を使います。
自分の実験用にしか使わないので、速度は遅くてもぜんぜん問題はありません。
やはり、有理数計算クラスを作るしかないみたいですね。
頑張って勉強してみます。
学歴っておいしい？さん、コードを書いていただいてありがとうございます。
参考にしてみます。

36 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/26(木) 19:33

>点列を色々な曲線で補間して、その滑らかさを比較する、
>ということをしています。点列を曲線で補間するときに
>連立方程式を使います。
ふむふむ。数学的な議論で終始するなら、連立方程式の厳密評価でも
イイかもしんないですね。

でももしアルゴリズムの研究なら「そのアルゴリズムがどのくらい
簡単に使えるか」も、アルゴリズムの評価の一端になると思います。
「多くの場合、非常に特異に近い係数行列になりますが、
その連立方程式が厳密に評価できれば良い補間です」とかだと、
実用上は使い物になりませんよね。
一般的な(陰的ピボット選択付きのLU分解などの)連立方程式
解法による評価も併せて行うべきだと思います。

それから、新規の手法の評価に有理数計算クラスをつくるなら、
徹底的にテストしないと、手法に不備があるのか、計算クラスに
不備があるのかで、わけわかんなくなるので要注意。

がんばってください。

37 名前： 名無しさん＠お腹いっぱい。 投稿日： 2000/10/26(木) 22:28

こっちのコードのほうが、参考になるだろう。。

http://www.interq.or.jp/jazz/iijima/equations/equations00.html