3 Windows フォーム・アプリケーション
4.5 連立一次方程式の解法
数学Cでは連立一次方程式を解くガウスの消去法が説明されBASIC のプログラムが与えられ ている。つぎのBASICのプログラムはガウスの消去法で連立一次方程式を解くものである。
100 REM ガウスの消去法 110 INPUT "N=";N 120 DIM A(N,N),B(N) 130 REM 方程式の入力 140 FOR I=1 TO N 150 FOR J=1 TO N
160 PRINT "A(";USING "#" I;USING "#" J;")=";:INPUT A(I,J) 170 NEXT J
180 PRINT "B(";USING "#" I;")=";:INPUT B(I) 190 NEXT I
200 REM 前進消去 210 FOR K=1 TO N-1 220 REM 行の入れ換え
230 AMAX=ABS(A(K,K)):KMAX=K 240 FOR I=K+1 TO N
250 IF ABS(A(I,K))>AMAX THEN KMAX=I:AMAX=ABS(A(I,K)) 260 NEXT I
270 IF AMAX<0.00001 THEN PRINT "解けない":END 280 IF KMAX=K THEN 340
290 FOR J=K TO N
300 T=A(K,J):A(K,J)=A(KMAX,J):A(KMAX,J)=T 310 NEXT J
320 T=B(K):B(K)=B(KMAX):B(KMAX)=T 330 REM 前進消去の中心部
340 FOR I=K+1 TO N
350 A(I,K)=A(I,K)/A(K,K) 360 FOR J=K+1 TO N
370 A(I,J)=a(I,J)-A(I,K)*A(K,J) 380 NEXT J
390 B(I)=B(I)-A(I,K)*B(K) 400 NEXT I
410 NEXT K
420 IF ABS(A(N,N))<0.00001 THEN PRINT "解けない" : END 430 REM 後退代入
440 FOR K=N TO 1 STEP -1 450 FOR I=K+1 TO N
460 B(K)=B(K)-A(K,I)*B(I) 470 NEXT
480 B(K)=B(K)/A(K,K) 490 NEXT K
500 REM 解の出力 510 PRINT "X ";
520 FOR K=1 TO N 530 PRINT B(K);
540 NEXT K 550 END
実行すると N=3
A( 1 1 )= ? 1 A( 1 2 )= ? 2 A( 1 3 )= ? -1 B( 1 )= ? -3
A( 2 1 )= ? 2 A( 2 2 )= ? 4 A( 2 3 )= ? 1 B( 2 )= ? 0 A( 3 1 )= ? 3 A( 3 2 )= ? 8 A( 3 3 )= ? 1 B( 3 )= ? -3
X 1.0000000000000007 -1.0000000000000002 2.0
となる。上の BASICのプログラムをC++に翻訳すると次のようになります。
#include <iostream>
#include <stdio.h>
#include <math.h>
using namespace std;
int main() { int n;
cout << "n="; cin >> n;
double **A, *B;
A = new double*[n];
for (int i=0; i<n; i++) A[i] = new double[n];
B = new double[n];
for (int i=0; i<n; i++) { for (int j=0; j<n; j++) {
cout << "A[" << i << "][" << j << "]="; cin >> A[i][j];
}
cout << "B[" << i << "]="; cin >> B[i];
}
for (int k=0; k<n; k++) {
double AMAX = fabs(A[k][k]); int KMAX = k;
for (int i=k+1; i<n; i++) { if (fabs(A[i][k]) > AMAX) {
KMAX = i;
AMAX = fabs(A[i][k]);
} }
if (AMAX < 0.00001) {
cout << "解けない" << endl;
for (int i=0; i<n; i++)
delete A[i];
delete A;
getchar();
getchar();
return 1;
}
if (KMAX != k) {
for (int j=k; j<n; j++) { double T = A[k][j];
A[k][j]= A[KMAX][j];
A[KMAX][j] = T;
}
double T = B[k];
B[k] = B[KMAX];
B[KMAX] = T;
}
for (int i=k+1; i<n; i++) { A[i][k] = A[i][k] / A[k][k];
for (int j=k+1; j<n; j++) {
A[i][j] = A[i][j] - A[i][k] * A[k][j];
}
B[i] = B[i] - A[i][k] * B[k];
} }
if (fabs(A[n-1][n-1]) < 0.00001) { cout << "解けない" << endl;
for (int i=0; i<n; i++) delete A[i];
delete A;
delete B;
getchar();
getchar();
return 1;
}
for (int k=n-1; k>=0; k--) { for (int i=k+1; i<n; i++) {
B[k] = B[k] - A[k][i] * B[i];
}
B[k] = B[k] / A[k][k];
}
cout << " X ";
for (int k=0; k<n; k++) { cout << B[k] << " ";
}
cout << endl;
for (int i=0; i<n; i++) delete A[i];
delete A;
delete B;
getchar();
getchar();
return 0;
}
ｎ次正方行列の逆行列を求めるには連立一次方程式をｎ個同時に解けばよい。次は逆行列を求め る BASICプログラムである。
100 REM ガウスの消去法 110 INPUT "N=";N 120 DIM A(N,N),B(N,N) 130 REM 行列の入力 140 FOR I=1 TO N 150 FOR J=1 TO N
160 PRINT "A(";USING "#" I;USING "#" J;")=";:INPUT A(I,J) 170 IF I=J THEN B(I,J)=1 ELSE B(I,J)=0
180 NEXT J 190 NEXT I 200 REM 前進消去 210 FOR K=1 TO N-1 220 REM 行の入れ換え
230 AMAX=ABS(A(K,K)):KMAX=K 240 FOR I=K+1 TO N
250 IF ABS(A(I,K))>AMAX THEN KMAX=I:AMAX=ABS(A(I,K)) 260 NEXT I
270 IF AMAX<0.00001 THEN PRINT "解けない":END 280 IF KMAX=K THEN 340
290 FOR J=K TO N
300 T=A(K,J):A(K,J)=A(KMAX,J):A(KMAX,J)=T 310 NEXT J
317 FOR J=1 TO N
320 T=B(K,J):B(K,J)=B(KMAX,J):B(KMAX,J)=T 323 NEXT J
330 REM 前進消去の中心部 340 FOR I=K+1 TO N
350 A(I,K)=A(I,K)/A(K,K) 360 FOR J=K+1 TO N
370 A(I,J)=a(I,J)-A(I,K)*A(K,J)
380 NEXT J 387 FOR J=1 TO N
390 B(I,J)=B(I,J)-A(I,K)*B(K,J) 393 NEXT J
400 NEXT I 410 NEXT K
420 IF ABS(A(N,N))<0.00001 THEN PRINT "解けない" : END 430 REM 後退代入
440 FOR K=N TO 1 STEP -1 450 FOR I=K+1 TO N 457 FOR J=1 TO N
460 B(K,J)=B(K,J)-A(K,I)*B(I,J) 463 NEXT J
470 NEXT
477 FOR J=1 TO N
480 B(K,J)=B(K,J)/A(K,K) 483 NEXT J
490 NEXT K 500 REM 解の出力 510 PRINT "逆行列 "
520 FOR K=1 TO N 527 FOR J=1 TO N 530 PRINT B(K,J);
533 NEXT J 536 PRINT 540 NEXT K 550 END
実行すると N=2
A( 1 1 )= ? 2 A( 1 2 )= ? 1 A( 2 1 )= ? 3 A( 2 2 )= ? 4 逆行列
0.8000000000000002 -0.20000000000000004 -0.6000000000000001 0.4
となる。
上のBASICのプログラムをC++に翻訳すると次のようになります。
#include <iostream>
#include <stdio.h>
#include <math.h>
using namespace std;
int main() { int n;
cout << "n="; cin >> n;
double **A, **B;
A = new double*[n];
for (int i=0; i<n; i++) A[i] = new double[n];
B = new double*[n];
for (int i=0; i<n; i++) B[i] = new double[n];
for (int i=0; i<n; i++) { for (int j=0; j<n; j++) {
cout << "A[" << i << "][" << j << "]="; cin >> A[i][j];
if (i == j) B[i][j] = 1;
else
B[i][j] = 0;
} }
for (int k=0; k<n; k++) {
double AMAX = fabs(A[k][k]); int KMAX = k;
for (int i=k+1; i<n; i++) { if (fabs(A[i][k]) > AMAX) {
KMAX = i;
AMAX = fabs(A[i][k]);
} }
if (AMAX < 0.00001) {
cout << "解けない" << endl;
for (int i=0; i<n; i++) delete A[i];
delete A;
for (int i=0; i<n; i++) delete B[i];
delete B;
getchar();
getchar();
return 1;
}
if (KMAX != k) {
for (int j=k; j<n; j++) { double T = A[k][j];
A[k][j]= A[KMAX][j];
A[KMAX][j] = T;
}
for (int j=k; j<n; j++) { double T = B[k][j];
B[k][j]= B[KMAX][j];
B[KMAX][j] = T;
} }
for (int i=k+1; i<n; i++) { A[i][k] = A[i][k] / A[k][k];
for (int j=k+1; j<n; j++) {
A[i][j] = A[i][j] - A[i][k] * A[k][j];
}
for (int j=k+1; j<n; j++) {
B[i][j] = B[i][j] - A[i][k] * B[k][j];
} } }
if (fabs(A[n-1][n-1]) < 0.00001) { cout << "解けない" << endl;
for (int i=0; i<n; i++) delete A[i];
delete A;
for (int i=0; i<n; i++) delete B[i];
delete B;
getchar();
getchar();
return 1;
}
for (int k=n-1; k>=0; k--) { for (int i=k+1; i<n; i++) {
for (int j=0; j<n; j++) {
B[k][j] = B[k][j] - A[k][i] * B[i][j];
} }
for (int j=0; j<n; j++) {
B[k][j] = B[k][j] / A[k][k];
}
}
cout << " 逆行列： " << endl;
for (int k=0; k<n; k++) { for (int j=0; j<n; j++) {
cout << B[k][j] << " ";
}
cout << endl;
}
for (int i=0; i<n; i++) delete A[i];
delete A;
for (int i=0; i<n; i++) delete B[i];
delete B;
getchar();
getchar();
return 0;
}
演習問題
1. 基本的な考え方はガウスの消去法と同じであるが、前進消去の段階で、対角線より下にある 成分だけでなく、対角線より上にある成分も同時に消去してしまう方法もある。これをガウ ス・ジョルダンの消去法あるいは掃出し法という。ガウス・ジョルダンの消去法により連立 一次方程式を解くプログラムを書け。
挑戦問題
1. バイオリズムは、身体の周期が２３日、感情の周期が２８日、知性の周期が３３日で、次の ように計算されます。
physical= sin(T /23∗2∗P I) sensitivity= sin(T /28∗2∗P I) intellectual= sin(T /33∗2∗P I)
ここで、Tは誕生日から計算する日までの日数で、PIは円周率です。各々の値が正ならば好 調、負ならば不調、零ならば好不調の変わり目の要注意日となります。ある日の前後のバイ オリズム曲線を描く場合には、上の式でTをその日の前後に変化させて、値を線で結べばよ いです。バイオリズムを表示するプログラムを作れ。
参考： sageによるプログラミング
sageという数式処理ソフトを使うとこのテキストで述べた計算が簡単にできます。
−x2−x3+x4= 0 x1+x2+x3+x4= 6 2x1+ 4x2+x3−2x4=−1 3x1+x2−2x3+ 2x4= 3 という連立方程式を解くには、sage のnotebookで
A = Matrix(QQ, [[0, -1, -1, 1],[1, 1, 1, 1],[2, 4, 1, -2][3, 1, -2, 2]]) B = vector([0, 6, -1, 3])
solution = A.solve_right(B) show(solution)
と入力し、
evaluateのボタンをクリックすれば良いです。
A=
2 5 4 3 1 2 5 4 6
(1)
の逆行列は
A = matrix(QQ, [[2, 5, 4],[3, 1, 2],[5, 4, 6]]) show(A.inverse())
と入力し、evaluateのボタンをクリックすれば良いです。
CやC++で不定積分を計算するのは大変難しいですが、sageでは簡単で var(’x’)
f(x) = e^x * cos(x)
f_int(x) = integrate(f, x) show(f_int)
と入力すればよく、√
(1−x2)の0 から1の定積分は f(x) = sqrt(1 - x^2)
f_integral = integrate(f, (x, 0, 1)) show(f_integral)
と入力すればいいです。