Vorheriger Artikel Portierung des Programms "FORTRAN77 Numerical Computing Programming" nach C und Python (Teil 2)
Ich habe den Code von Kapitel 1 von Teil 1, Teil 2 und Fortran 77 Numerical Calculation Programming überprüft. In dem Buch gibt es danach drei Programme in Kapitel 1, aber alle prüfen, was der Fehler sein wird. Wie mir @cure_honey letztes Mal sagte, stimmt das Berechnungsergebnis unter Verwendung der Gleitkommazahl basierend auf der Hexadezimalzahl nicht mit dem Ergebnis mit der Gleitkommazahl der aktuellen IEEE754-Methode überein und weicht von der Überprüfung im Buch ab. Ich werde den Code hier einfach überspringen. Wenn Sie den Fehler überspringen, folgt die Geschichte der Zufallszahlen. Da es sowohl in C als auch in Python Funktionen zum Generieren von Zufallszahlen gibt, ist das nächste Kapitel 5, "Simultane lineare Gleichungen und LU-Zerlegung".
Die LU-Zerlegung wird durch Zeichnen der Formeln aus Kapitel 5, S. 52 des Buches zusammengefasst.
$ n \ times 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}
Und wenn $ b $ und $ x $ die folgenden Vektoren sind,
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}
\\
Dafür $ A $ und $ b $
Ax = b\tag{4}
Betrachten Sie ein lineares Gleichungssystem für $ x $.
Die LU-Zerlegung ist das Produkt der unteren linken Dreiecksmatrix $ L $ und der oberen rechten Dreiecksmatrix $ U $ für eine gegebene Matrix $ A $. Mit anderen Worten
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
Es ist zu bestimmen, ■ dass es wird. In Wirklichkeit jedoch
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}
Es ist auf die Form beschränkt. Auf diese Weise wird die Gleichung $ (4) $
LUx=b\tag{6}
Und um dies zu lösen, zuerst
Ly=b\tag{7}
Finden Sie die Lösung $ y $ von und dann die Gleichung mit dieser Lösung $ y $ auf der rechten Seite
Ux=y\tag{8}
Alles was Sie tun müssen, ist zu lösen, um $ x $ zu finden. Es scheint ärgerlich zu sein, weil die Anzahl der Gleichungen auf zwei gestiegen ist, aber wenn die Matrix zu einer dreieckigen Matrix wird, ist es viel einfacher zu lösen.
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 Hauptprogramm. Die Anweisung PARAMETER definiert NDIM, das die Größe des Arrays für das Beispiel angibt.
ludata Zunächst das Unterprogramm ludata
N=10
Wann
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}
Ich mache eine Linie namens. Zu diesem Zeitpunkt ist die Matrix $ b $
b_i = \sum_{j=1}^na_{ij}j
Wird gesucht. Was ich tun möchte, ist für diese $ A $ und $ b $
Ax = b
War $ x $ zu finden.
ancomp Als nächstes wird das Unterprogramm ancomp aufgerufen
||A||_1 = \max_{1\leq k\leq n} \sum_{i=1}^n|a_{ik}|\tag{9}
Die Anzahl der Bedingungen wird basierend auf der Norm der Matrix festgelegt. Dies wird in einer Variablen namens $ ANORM $ festgelegt. Der Grund, warum die Arbeitsvariable $ S $ nicht als reeller Zahlentyp mit doppelter Genauigkeit deklariert wird, ist, dass diese Norm selbst nicht viel System erfordert.
decomp Eine Unterroutine, in der die Dekomposition die LU-Zerlegung des Hauptthemas durchführt. Die Argumente sind das zweidimensionale Array $ A $, das der Matrix $ A $ entspricht, das eindimensionale ganzzahlige Array $ IP $, das dem Indexvektor $ p $ entspricht, und die Matrixdimension $ N $. Zunächst das Array $ IP $,
IP(k) = k, \quad k=1,2,\cdots, n
Führen Sie nach dem Initialisieren mit eine LU-Zerlegung nach der Gaußschen Eliminierungsmethode durch, während Sie eine teilweise Auswahl des Pivots durchführen.
solve Solve ist eine Unterroutine, die die Gleichung $ (4) $ unter Verwendung des Ergebnisses der LU-Zerlegung löst.
estcnd Estcnd ist eine Unterroutine, die eine Schätzung der Anzahl der Bedingungen in der Matrix $ A $ berechnet. Es wird jedoch angenommen, dass die Matrix $ A $ im Voraus durch Zerlegung LU-zerlegt wurde.
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;
}
Ersetzen Sie einfach ehrlich den Fortran-Code. Ich fragte mich, ob ich die Variablen global definieren könnte, aber ich persönlich konnte es nicht tun.
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 hat eine leistungsstarke numerische Berechnungsbibliothek namens SciPy, und natürlich gibt es auch eine LU-Zerlegungsfunktion ... Deshalb habe ich nur Daten geschrieben, als ich Daten gemacht habe.
Wenn Sie danach suchen, frage ich mich, ob es auch eine schöne Bibliothek in C-Sprache gibt. Vorerst war es ziemlich ärgerlich, die Formel zum dritten Mal zu schreiben, als ich sie lange Zeit niederlegte.