Article précédent Portage du programme "FORTRAN77 Numerical Computing Programming" vers C et Python (partie 2)
J'ai vérifié le code du chapitre 1 de la partie 1, partie 2 et Programmation de calcul numérique Fortran 77. Dans le livre, il y a trois programmes après cela dans le chapitre 1, mais tous vérifient quelle sera l'erreur. Comme @cure_honey me l'a dit la dernière fois, le résultat du calcul utilisant le nombre à virgule flottante basé sur le nombre hexadécimal ne correspond pas au résultat avec le nombre à virgule flottante de la méthode IEEE754 actuelle, et il s'écarte de la vérification dans le livre. Je vais juste sauter le code ici. Si vous ignorez l'erreur, voici l'histoire des nombres aléatoires. Puisqu'il existe des fonctions à la fois en C et en Python pour générer des nombres aléatoires, si vous l'ignorez également, voici le chapitre 5, "Equations linéaires simultanées et décomposition LU".
La décomposition LU est résumée en tirant les formules du chapitre 5, P.52 du livre.
$ n \ fois n $ Matrix Matrix $ A $
A=\begin{pmatrix}
a_{11} & a_{12} & \cdots & a_{1n} \\
a_{21} & a_{22} & \cdots & a_{2n} \\
\vdots & \vdots & & \vdots \\
a_{n1} & a_{n2} & \cdots & a_{nn}
\end{pmatrix}\tag{1}
Et lorsque $ b $ et $ x $ sont respectivement les vecteurs suivants,
b = \begin{pmatrix}
b_1\\
b_2\\
\vdots\\
b_n
\end{pmatrix}\tag{2}
\\
x = \begin{pmatrix}
x_1\\
x_2\\
\vdots\\
x_n
\end{pmatrix}\tag{3}
\\
Pour ces $ A $ et $ b $
Ax = b\tag{4}
Considérons un système d'équations linéaires pour $ x $.
La décomposition LU est le produit de la matrice triangulaire inférieure gauche $ L $ et de la matrice triangulaire supérieure droite $ U $ pour une matrice donnée $ A $. En d'autres termes
A=\begin{pmatrix}
■&0&0&0&0\\
■&■&0&0&0\\
■&■&■&0&0\\
■&■&■&■&0\\
■&■&■&■&■
\end{pmatrix}
\begin{pmatrix}
■&■&■&■&■\\
0&■&■&■&■\\
0&0&■&■&■\\
0&0&0&■&■\\
0&0&0&0&■
\end{pmatrix}
\equiv
LU
Il s'agit de déterminer ■ pour qu'il devienne. Cependant, en réalité
A=\begin{pmatrix}
1&0&0&0&0\\
■&1&0&0&0\\
■&■&1&0&0\\
■&■&■&1&0\\
■&■&■&■&1
\end{pmatrix}
\begin{pmatrix}
■&■&■&■&■\\
0&■&■&■&■\\
0&0&■&■&■\\
0&0&0&■&■\\
0&0&0&0&■
\end{pmatrix}
\equiv
LU\tag{5}
Il est limité à la forme. En faisant cela, l'équation $ (4) $
LUx=b\tag{6}
Pour résoudre ce problème, d'abord
Ly=b\tag{7}
Trouvez la solution $ y $ de, puis l'équation avec cette solution $ y $ sur le côté droit
Ux=y\tag{8}
Tout ce que vous avez à faire est de résoudre pour trouver $ x $. Cela semble gênant car le nombre d'équations est passé à deux, mais lorsque la matrice devient une matrice triangulaire, c'est beaucoup plus facile à résoudre.
sdecom.f
PROGRAM SDECOM
* *
* Sample program of DECOMP and SOLVE *
* *
PARAMETER (NDIM = 20)
*
REAL A(NDIM,NDIM),B(NDIM),X(NDIM)
REAL W(NDIM),V(NDIM)
INTEGER IP(NDIM)
*
DO 10 NN = 1, 2
*
N = 10 * NN
*
CALL LUDATA (NDIM, N, A, B)
*
CALL ANCOMP (NDIM, N, A, ANORM)
*
WRITE (*,2000) ANORM
2000 FORMAT (/' ---- DECOMP and SOLVE ----'/
$ ' 1-norm of A = ',1PE13.6)
CALL DECOMP (NDIM, N, A, W, IP)
CALL SOLVE (NDIM, N, A, B, X, IP)
CALL ESTCND (NDIM, N, A, V, IP, W,
$ ANORM, COND)
*
WRITE (*,2001) (X(I), I=1,N)
2001 FORMAT (' solution X(I)'/(1X,5F9.5))
*
WRITE (*,2002) COND
2002 FORMAT (' condition number = ',1PE13.6)
*
10 CONTINUE
*
END
ludata.f
SUBROUTINE LUDATA (NDIM, N, A, B)
* *
* Sample data of DECOMP and SOLVE *
* *
REAL A(NDIM,N),B(N)
DO 510 I =1, N
DO 520 J = 1, N
A(I,J) = 0.0
520 CONTINUE
510 CONTINUE
*
A(1,1) = 4.0
A(1,2) = 1.0
DO 530 I = 2, N - 1
A(I,I-1) = 1.0
A(I,I) = 4.0
A(I,I+1) = 1.0
530 CONTINUE
A(N,N-1) = 1.0
A(N,N) = 4.0
DO 540 I = 1, N
B(I) = 0.0
DO 550 J = 1, N
B(I) = B(I) + A(I,J) * J
550 CONTINUE
540 CONTINUE
*
RETURN
END
ancomp.f
SUBROUTINE ANCOMP (NDIM, N, A, ANORM)
* *
* Compute 1-norm of A toestimate *
* condition number of A *
* *
REAL A(NDIM,N)
*
ANORM = 0.0
DO 10 K = 1, N
S = 0.0
DO 520 I = 1, N
S = S + ABS(A(I,K))
520 CONTINUE
IF (S .GT. ANORM) ANORM = S
10 CONTINUE
*
RETURN
END
decomp.f
SUBROUTINE DECOMP (NDIM, N, A, W, IP)
* *
* LU decomposition of N * N matrix A *
* *
REAL A(NDIM,N),W(N)
INTEGER IP(N)
*
DATA EPS / 1.0E-75 /
*
DO 510 K = 1, N
IP(K) = K
510 CONTINUE
*
DO 10 K = 1, N
L = K
AL = ABS(A(IP(L),K))
DO 520 I = K+1, N
IF (ABS(A(IP(I),K)) .GT. AL) THEN
L = I
AL = ABS(A(IP(L),K))
END IF
520 CONTINUE
*
IF (L .NE. K) THEN
*
* ---- row exchange ----
*
LV = IP(K)
IP(K) = IP(L)
IP(L) = LV
END IF
*
IF (ABS(A(IP(K),K)) .LE. EPS) GO TO 900
*
* ---- Gauss elimination ----
*
A(IP(K),K) = 1.0 / A(IP(K),K)
DO 30 I = K+1, N
A(IP(I),K) = A(IP(I),K) * A(IP(K),K)
DO 540 J = K+1, N
W(J) = A(IP(I),J) - A(IP(I),K) * A(IP(K),J)
540 CONTINUE
DO 550 J = K+1, N
A(IP(I),J) = W(J)
550 CONTINUE
30 CONTINUE
*
10 CONTINUE
*
RETURN
*
* ---- matrix singular ----
*
900 CONTINUE
WRITE (*,2001) K
2001 FORMAT (' (DECOMP) matrix singular ( at ',I4,
$ '-th pivot)')
N = - K
RETURN
*
END
solve.f
SUBROUTINE SOLVE (NDIM, N, A, B, X, IP)
* *
* Solve system of linear equations *
* A * X = B *
* *
REAL A(NDIM,N),B(N),X(N)
INTEGER IP(N)
*
REAL*8 T
*
* ---- forward substitution ----
*
DO 10 I = 1,N
T = B(IP(I))
DO 520 J = 1, I-1
T = T - A(IP(I),J) * DBLE(X(J))
520 CONTINUE
X(I) = T
10 CONTINUE
*
* ---- backward substitution ----
*
DO 30 I = N, 1, -1
T = X(I)
DO 540 J = I+1, N
T = T - A(IP(I),J) * DBLE(X(J))
540 CONTINUE
X(I) = T * A(IP(I),I)
30 CONTINUE
*
RETURN
END
estcnd.f
SUBROUTINE ESTCND(NDIM, N, A, V, IP, Y,
$ ANORM, COND)
* *
* Estimate condition number of A *
* *
REAL A(NDIM,N),V(N),Y(N)
INTEGER IP(N)
*
REAL*8 T
*
DO 10 K = 1,N
T = 0.0
DO 520 I = 1, k-1
T = T - A(IP(I),K) * DBLE(V(I))
520 CONTINUE
S = 1.0
IF (T .LT. 0.0) S = -1.0
V(K) = (S + T) * A(IP(K),K)
10 CONTINUE
*
DO 30 K = N, 1, -1
T = V(K)
DO 530 I = K+1, N
T = T - A(IP(I),K) * DBLE(Y(IP(I)))
530 CONTINUE
Y(IP(K)) = T
30 CONTINUE
*
YMAX = 0.0
DO 540 I = 1, N
IF (ABS(Y(I)) .GT. YMAX) YMAX = ABS(Y(I))
540 CONTINUE
*
COND = ANORM * YMAX
*
RETURN
END
sdecom Programme principal. L'instruction PARAMETER définit NDIM, qui donne la taille du tableau pour l'exemple.
ludata Tout d'abord, le sous-programme appelé ludata
N=10
Quand
A=\begin{pmatrix}
4.0 & 1.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\
1.0 & 4.0 & 1.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\
0.0 & 1.0 & 4.0 & 1.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\
0.0 & 0.0 & 1.0 & 4.0 & 1.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 \\
0.0 & 0.0 & 0.0 & 1.0 & 4.0 & 1.0 & 0.0 & 0.0 & 0.0 & 0.0 \\
0.0 & 0.0 & 0.0 & 0.0 & 1.0 & 4.0 & 1.0 & 0.0 & 0.0 & 0.0 \\
0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1.0 & 4.0 & 1.0 & 0.0 & 0.0 \\
0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1.0 & 4.0 & 1.0 & 0.0 \\
0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 4.0 & 1.0 \\
0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 1.0 & 4.0 \\
\end{pmatrix}
b=\begin{pmatrix}
6.00000 & 12.00000 & 18.00000 & 24.00000 & 30.00000 & 36.00000 & 42.00000 & 48.00000 & 54.00000 & 49.00000
\end{pmatrix}
Je fais une ligne appelée. A ce moment, la matrice $ b $ est
b_i = \sum_{j=1}^na_{ij}j
Est recherché. Ce que je veux faire, c'est pour ces $ A $ et $ b $
Ax = b
Était de trouver $ x $.
ancomp Ensuite, le sous-programme appelé ancomp
||A||_1 = \max_{1\leq k\leq n} \sum_{i=1}^n|a_{ik}|\tag{9}
Le nombre de conditions est fixé en fonction de la norme de la matrice. Ceci est défini dans une variable appelée $ ANORM $. La raison pour laquelle la variable de travail $ S $ n'est pas déclarée ici comme un type de nombre réel à double précision est que cette norme elle-même ne nécessite pas beaucoup de système.
decomp Un sous-programme dans lequel decomp effectue la décomposition LU du sujet principal. Les arguments sont le tableau bidimensionnel $ A $ correspondant à la matrice $ A $, le tableau d'entiers unidimensionnel $ IP $ correspondant au vecteur d'index $ p $ et la dimension matricielle $ N $. Tout d'abord, le tableau $ IP $,
IP(k) = k, \quad k=1,2,\cdots, n
Après l'initialisation avec, effectuez la décomposition LU par la méthode d'élimination gaussienne tout en effectuant une sélection partielle du pivot.
solve Solve est un sous-programme qui résout l'équation $ (4) $ en utilisant le résultat de la décomposition LU.
estcnd Estcnd est un sous-programme qui calcule une estimation du nombre de conditions dans la matrice $ A $. Cependant, on suppose que la matrice $ A $ a été décomposée par LU par décompression à l'avance.
sdecom.h
#ifndef NDIM
#define NDIM 20
#endif
void ludata(int n, float a[NDIM][NDIM], float b[NDIM]);
float ancomp(int n, float a[NDIM][NDIM]);
int decomp(int n, float a[NDIM][NDIM], float w[NDIM], int ip[NDIM]);
void solve(int n, float a[NDIM][NDIM], float b[NDIM], float x[NDIM], int ip[NDIM]);
float estcnd(int n, float a[NDIM][NDIM], float v[NDIM], int ip[NDIM], float y[NDIM], float anorm);
sdecom.c
#include <stdio.h>
#include<stdlib.h>
#include "sdecom.h"
int main(void)
{
float a[NDIM][NDIM], b[NDIM], x[NDIM];
float w[NDIM], v[NDIM];
int ip[NDIM];
int nn;
int n;
float anorm, cond;
int i,j;
for(nn = 1; nn < 3; nn++){
n = 10 * nn;
ludata(n, a, b);
anorm = ancomp(n, a);
printf(" ---- DECOMP and SOLVE ----\n 1-norm of A = %13.6E\n"
, anorm);
n = decomp(n, a, w, ip);
solve(n, a, b, x, ip);
cond = estcnd(n, a, v, ip, w, anorm);
printf(" solution X(I)\n");
for(i=0; i<n; i++){
if(i != 0 && i % 5 == 0) {
printf("\n");
}
printf(" %.5f", x[i]);
}
printf("\n condition number = %13.6E\n", cond);
}
return 0;
}
ludata.c
#include <stdio.h>
#include "sdecom.h"
void ludata(int n, float a[NDIM][NDIM], float b[NDIM])
{
int i, j;
for(i=0; i<n; i++){
for(j=0; j<n; j++){
a[i][j] = 0.0f;
}
}
a[0][0] = 4.0f;
a[0][1] = 1.0f;
for(i=1; i<n-1; i++){
a[i][i-1] = 1.0f;
a[i][i] = 4.0f;
a[i][i+1] = 1.0f;
}
a[n-1][n-2] = 1.0f;
a[n-1][n-1] = 4.0f;
for(i=0; i<n; i++){
b[i] = 0.0f;
for(j=0; j<n; j++) {
b[i] += (a[i][j] * (j+1));
}
}
}
ancomp.c
#include <stdio.h>
#include <math.h>
#include "sdecom.h"
float ancomp(int n, float a[NDIM][NDIM])
{
int i,k;
float s;
float anorm = 0.0f;
for(k=0; k<n; k++){
s = 0.0f;
for(i=0; i<n; i++) {
s = s + fabsf(a[i][k]);
}
if (s > anorm){
anorm = s;
}
}
return anorm;
}
decomp.c
#include <stdio.h>
#include <math.h>
#include <float.h>
#include "sdecom.h"
int decomp(int n, float a[NDIM][NDIM], float w[NDIM], int ip[NDIM])
{
//double eps = 1.0e-75;
int i,j,k,l;
float al;
int lv;
for(k=0; k<n; k++){
ip[k] = k;
}
for(k=0; k<n; k++){
l = k;
al = fabsf(a[ip[l]][k]);
for(i=k+1; i<n; i++){
if(fabsf(a[ip[i]][k]) > al) {
l = i;
al = fabsf(a[ip[l]][k]);
}
}
if(l != k) {
lv = ip[k];
ip[k] = ip[l];
ip[l] = lv;
}
if((fabsf(a[ip[k]][k]) <= DBL_EPSILON)) {
printf(" (DECOMP) matrix singular ( at %4d-th pivot)\n", k);
n = -k;
break;
} else {
a[ip[k]][k] = 1.0f / a[ip[k]][k];
for(i=k+1; i<n; i++) {
a[ip[i]][k] = a[ip[i]][k] * a[ip[k]][k];
for(j=k+1; j<n; j++) {
w[j] = a[ip[i]][j] - a[ip[i]][k] * a[ip[k]][j];
}
for(j=k+1; j<n; j++) {
a[ip[i]][j] = w[j];
}
}
}
}
return n;
}
solve.c
#include <stdio.h>
#include "sdecom.h"
void solve(int n, float a[NDIM][NDIM], float b[NDIM], float x[NDIM], int ip[NDIM])
{
double t;
int i,j;
for(i=0; i<n; i++){
x[i] = 0.0f;
}
for(i=0; i<n; i++){
t = b[ip[i]];
for(j=0; j<i; j++){
t = t - a[ip[i]][j] * (double)x[j];
}
x[i] = t;
}
for(i=n-1; i>-1; i--){
t = x[i];
for(j=i+1; j <n; j++){
t = t - a[ip[i]][j] * (double)x[j];
}
x[i] = t * a[ip[i]][i];
}
}
estcnd.c
#include <stdio.h>
#include <math.h>
#include "sdecom.h"
float estcnd(int n, float a[NDIM][NDIM], float v[NDIM], int ip[NDIM], float y[NDIM], float anorm)
{
float cond;
double t;
int i,k;
float s;
float ymax;
for(k=0; k<n; k++){
t = 0.0;
for(i=0; i<k; i++){
t = t - a[ip[i]][k] * (double)v[i];
}
s = 1.0f;
if( t < 0.0f) {
s = -1.0f;
}
v[k] = (s + t) * a[ip[k]][k];
}
for(k=n-1; k>-1; k--){
t = v[k];
for(i=k+1; i<n; i++){
t = t - a[ip[i]][k] * (double)y[ip[i]];
}
y[ip[k]] = t;
}
ymax = 0.0f;
for(i=0; i<n; i++){
if(fabsf(y[i]) > ymax){
ymax = fabsf(y[i]);
}
}
cond = anorm * ymax;
return cond;
}
Remplacez simplement le code Fortran. Je me suis demandé si je pouvais définir les variables globalement, mais personnellement, je ne pouvais pas le faire.
import scipy.linalg as linalg
import numpy as np
def ludata(n):
a = [[0.0] * n for i in range(n) ]
a[0][0] = 4.0
a[0][1] = 1.0
for i in range(1,n-1):
a[i][i-1] = 1.0
a[i][i] = 4.0
a[i][i+1] = 1.0
a[n-1][n-2] = 1.0
a[n-1][n-1] = 4.0
b = [0.0] * n
for i in range(n):
for j in range(n):
b[i] += (a[i][j] * (j+1))
return a,b
for nn in range(1,3):
n = 10 * nn
a,b = ludata(n)
LU = linalg.lu_factor(a)
x = linalg.lu_solve(LU, b)
print(x)
Python dispose d'une puissante bibliothèque de calcul numérique appelée SciPy, et bien sûr, il existe également une fonction de décomposition LU ... C'est pourquoi je n'écrivais des données que lorsque j'en faisais.
Si vous le cherchez, je me demande s'il existe également une belle bibliothèque publiée en langage C. Pour le moment, c'était assez ennuyeux d'écrire la formule pour la troisième fois que je l'ai posée pendant longtemps.
Recommended Posts