「ガウスの消去法」というのがありますが、ゼロ除算が怖かったです。
まず、自分で係数を文字で置いて実際に解いてみました。
この方法で公式を作っちゃえばいいんじゃないか?と思ったのです。
しかし、大変でした。
そこで青チャート(数学IIIC)を眺めていたら・・・
こんなものがあったのです。
(ちゃららら)クラーメルの公式~
これを使うとね、簡単に連立方程式の解が求まるんだよ。(CV:大山のぶ代)←またお前か
二元連立方程式と三元連立方程式の例がありましたが、もしかしてこれはn元連立方程式に応用できるのでは?
ということでプログラムを組んでみました。
n×n行列(n≧4)の行列式の求め方は書いてありませんでしたが、家にあった電子辞書で調べたら解決しました。
まずは検証用にn元連立方程式を作るプログラム。
nを入力すると自動で生成してくれます。
#include
#include
#include
int main(void) {
int n;
int* ans;
int keisuu;
int sum;
int i,j;
scanf("%d",&n);
ans=malloc(sizeof(int)*n);
if(ans==NULL)return 1;
srand((unsigned int)time(NULL));
printf("%d\n",n);
for(i=0;i
#include
double** allocMatrix(int n) {
double** result;
int i;
/*列のメモリを確保*/
result=malloc(sizeof(double*)*n);
/*メモリ確保エラー*/
if(result==NULL)return NULL;
/*各行のメモリを確保*/
for(i=0;i=0;i--)free(result[i]);
free(result);
return NULL;
}
}
return result;
}
void freeMatrix(double** mat,int n) {
int i;
for(i=0;i=n)return 0;
/*置換する*/
if(j!=i) {
temp=tikan2[i];
tikan2[i]=tikan2[j];
tikan2[j]=temp;
tikannum++;
}
}
free(tikan2);
return tikannum%2;
}
double getOneMul(double** mat,int* tikan,int n) {
double result=1;
int i;
for(i=0;i=max)return getOneMul(mat,tikan,max);
/*全ての数字を試す*/
for(i=0;i=n) {
/*入れる*/
tikan[n]=i;
/*次の数字を試す*/
result+=getTikan(mat,tikan,n+1,max);
}
}
return result;
}
double getGyouretusiki(double** mat,int n) {
double result;
int* tikan;
/*メモリの確保*/
tikan=malloc(sizeof(int)*n);
if(tikan==NULL)return 0;
/*計算*/
result=getTikan(mat,tikan,0,n);
/*メモリの開放*/
free(tikan);
/*結果を返す*/
return result;
}
/*クラーメルの公式的なサムシング*/
int solveHouteisiki(double* result,double** mat,double* vec,int n) {
double** mat2;
int d;
int i,j,k;
/*そのまま行列式を求める*/
d=getGyouretusiki(mat,n);
/*エラー*/
if(d==0)return 0;
/*第二の行列を確保*/
mat2=allocMatrix(n);
/*エラー*/
if(mat2==NULL)return 0;
/*計算*/
for(i=0;i<n;i++) {
/*行列を作成*/
for(j=0;j<n;j++) {
for(k=0;k<n;k++) {
if(k==i) {
mat2[j][k]=vec[j];
} else {
mat2[j][k]=mat[j][k];
}
}
}
/*計算*/
result[i]=getGyouretusiki(mat2,n)/d;
}
/*メモリの開放*/
freeMatrix(mat2,n);
}
/*
使い方
まず、n元連立方程式を解く場合のnを整数で入力します。
次に、例えばn=3のとき、
一つの式がax+by+cz=dだとしたら、
a b c d
と入力します。
これをn回繰り返します。
すると、結果が最初の変数から順に表示されます。
この例だと
x
y
z
と表示されます。
入力のa、b、c、dと出力のx、y、zには具体的な実数が入ります。
*/
int main(void) {
int n;
double** mat;
double* vec;
double* result;
int i,j;
/*nの入力*/
scanf("%d",&n);
if(n<=0)return 1;
/*行列の確保*/
mat=allocMatrix(n);
if(mat==NULL)return 1;
/*ベクトルの確保*/
vec=allocVector(n);
if(vec==NULL) {
freeMatrix(mat,n);
return 1;
}
result=allocVector(n);
if(result==NULL) {
freeMatrix(mat,n);
freeVector(vec);
return 1;
}
/*値の入力*/
for(i=0;i<n;i++) {
for(j=0;j<n;j++) {
scanf("%lf",&mat[i][j]);
}
scanf("%lf",&vec[i]);
}
/*計算*/
solveHouteisiki(result,mat,vec,n);
/*結果の表示*/
for(i=0;i<n;i++) {
printf("%g\n",result[i]);
}
/*メモリの開放*/
freeMatrix(mat,n);
freeVector(vec);
freeVector(result);
/*正常終了*/
return 0;
}
1≦n≦8の時は(多分)成功!
n=8にすると少し重いですが、この範囲では正確に求まっています。
しかし、n=9にするとかなり重くなり、結果もめちゃくちゃです。(たまに正確、オーバーフロー?)
n=10の時はさらに重く、結果もめちゃくちゃです。
まあ、普通のプログラムで十元連立方程式なんて多分解かないだろうから、実用に耐える範囲でしょう。