XProfan

Español

p.specht

Der Householder-Algorithmus
=====================
En kritischen Datenkonstellationen versagen übliche Lösungsverfahren para Gleichungssysteme. Das Gauss'sche Eliminationsverfahren etwa produziert en großen Werte-Spannweiten innerhalb el Matrizenspalten erhebliche Fehler, el a absoluter Unbrauchbarkeit el Ergebnisse führen puede, sin dass el Anwender lo igual merken sería. Für solche Fälle hatte el US-Amerikanische Mathematiker Alston Householder 1958 el Concepto, el No-Diagonalwerte por geschickt gewählte Spiegelungen a Null a hacer. Besonders effizient gelingt en el Diagonalsymmetrischen Matrizen, como ellos en el Physik oft vorkommen. (Ansonsten puede ser jede beliebige Matrix por Multiplikation con ihrer Transponierten a uno Diagonalsymmetrischen Matrix hacer y el Resultat entonces LU-Zerlegen).

In dieser Demo para private Bildungszwecke wurde Fortran's Arrayindex-Base '1' beibehalten! In XProfan una ziemliche Stolperfalle!

Kompilieren Marca Separación

Título de la ventana "Inversion uno rellen, diagonalsymetrischen Matrix mittels Householder-Algorithmus"
'Solución por multidimensionaler Spiegelung, Algorithmus de Alston Scott Householder, USA 1958
'Orig.F90-Source: https://jean-pierre.moreau.pagesperso-orange.fr/Fortran/householder_f90.txt
'(D)Demo-Migration después de XProfan11.2a 2014-10 by P.Pájaro carpintero, Wien;  OHNE JEDWEDE GEWÄHR!
'Verfahrensdetails siehe https://de.wikipedia.org/wiki/Householdertransformation
'***************************************************************
'{  Inversion of a real square matrix by Householder's method  *
'
'* SAMPLE RUN [en el Original sólo single precision, Anm. P.Pájaro carpintero]*
'* 4                                                           *
'* 1.0       0.5        0.333333  0.25                         *
'* 0.5       0.333333   0.25      0.2                          *
'* 0.333333  0.25       0.2       0.166667                     *
'* 0.25      0.2        0.166667  0.142857                     *
'
'* Output data:
'* Entrada square matrix
'*  0.100000E+01   0.500000E+00   0.333333E+00   0.250000E+00  *
'*  0.500000E+00   0.333333E+00   0.250000E+00   0.200000E+00  *
'*  0.333333E+00   0.250000E+00   0.200000E+00   0.166667E+00  *
'*  0.250000E+00   0.200000E+00   0.166667E+00   0.142857E+00  *
'*                                                             *
'* Transformation # 1                                          *
'* -0.119315E+01  -0.670493E+00  -0.474932E+00  -0.369836E+00  *
'*  0.000000E+00   0.664812E-01   0.657297E-01   0.586884E-01  *
'*  0.000000E+00   0.720990E-01   0.771533E-01   0.724594E-01  *
'*  0.000000E+00   0.665741E-01   0.745319E-01   0.722012E-01  *
'* Transformation # 2                                          *
'* -0.119315E+01  -0.670493E+00  -0.474932E+00  -0.369836E+00  *
'*  0.000000E+00  -0.118533E+00  -0.125656E+00  -0.117542E+00  *
'*  0.000000E+00   0.000000E+00   0.257162E-02   0.378339E-02  *
'*  0.000000E+00   0.000000E+00   0.566533E-02   0.878781E-02  *
'* Transformation # 3                                          *
'* -0.119315E+01  -0.670493E+00  -0.474932E+00  -0.369836E+00  *
'*  0.000000E+00  -0.118533E+00  -0.125656E+00  -0.117542E+00  *
'*  0.000000E+00   0.000000E+00  -0.622167E-02  -0.956580E-02  *
'*  0.000000E+00   0.000000E+00   0.000000E+00   0.187201E-03  *
'
'* Inverted matrix:
'*  0.160327E+02  -0.120367E+03   0.240887E+03  -0.140579E+03  *
'* -0.120367E+03   0.120414E+04  -0.271001E+04   0.168653E+04  *
'*  0.240887E+03  -0.271001E+04   0.650422E+04  -0.421582E+04  *
'* -0.140579E+03   0.168653E+04  -0.421582E+04   0.281033E+04  *
'
'* Determinante: -0.164722E-06
'
'* AP=A*A^-1 matrix:                                           *
'*  0.100000E+01   0.341061E-12  -0.682121E-12   0.341061E-12  *
'*  0.710543E-14   0.100000E+01   0.000000E+00   0.000000E+00  *
'* -0.710543E-14   0.113687E-12   0.100000E+01   0.113687E-12  *
'* -0.355271E-14   0.284217E-13   0.000000E+00   0.100000E+01  *
'* ----------------------------------------------------------- *
'* Ref.: "Méthodes de calcul numérique - Tome 2 - By Claude    *
'*        Nowakowski, PSI Editions, 1984".                     *
'*                  Fortran 90 Versión By J-P Moreau, Paris.   *
'*                            (www.jpmoreau.fr)                *
'*       XProfan11.2a Migration by P.Pájaro carpintero, Vienna            *
'}
Ventana de Estilo 24
Ventana 0,0-%maxx,%maxy-40
Font 2
set("decimals",17)
Declarar A![25,50]
Imprimir "\n\n               DER HOUSHOLDER-ALGORITHMUS\n               =========================="
Householder(0,A![])' En Angabe n = 0 Dimensionen se el Testbeispiel berechnet
Waitinput
FIN

Proc Householder :parámetros n&,A![]

    declarar A0![25,25], V![25],D!,B![25],AI![25,25],AP![25,25],X![25],i&,j&,k&

    if n&=0'Read matrix to be inverted

        n&=4
        A![1,1]=1.0 :A![1,2]=0.5   :A![1,3]=0.333333333333333333:A![1,4]=0.25
        A![2,1]=0.5 :A![2,2]=0.333333333333333333 :A![2,3]=0.25 :A![2,4]=0.2
        A![3,1]=0.333333333333333333:A![3,2]=0.25 :A![3,3]=0.2  :A![3,4]=0.166666666666666666
        A![4,1]=0.25 :A![4,2]=0.2  :A![4,3]=0.166666666666666666:A![4,4]=0.142857142857142857

    endif

    Imprimir "\n es para diagonalsymetrische Ausgangsmatrizen vorgesehen,   BEISPIEL:\n"

    whileloop 1,n&:i&=&Loop

        whileloop 1,n&:j&=&Loop

            If j& = n&

                imprimir tab(j&*27);if(A![i&,j&]<0,""," ");format$("%g",A![i&,j&])

            Más

                imprimir tab(j&*27);if(A![i&,j&]<0,""," ");format$("%g",A![i&,j&]),

            endif

        endwhile

    endwhile

    whileloop 1,n&:i&=&Loop

        whileloop 1,n&:j&=&Loop

            A![i&,j&+n&] = 0
            A0![i&,j&] = A![i&,j&]

        endwhile

        A![i&,i&+n&] = 1

    endwhile

    'Transform A into triangular matrix
    S1000(n&,A![],V![])
    'N linear systems to solve

    whileloop 1,n&:k&=&Loop

        whileloop 1,n&:i&=&Loop

            B![i&] = A![i&,K&+n&]

        endwhile

        'Solve triangular system
        S2000(n&,A![],B![],X![])

        whileloop 1,n&:i&=&Loop

            AI![i&,K&]=X![i&]

        endwhile

    endwhile'k bucle

    imprimir "\n\n Invertierte Matrix:",
    ' Inverted matrix:'

    whileloop 1,n&:i&=&Loop

        whileloop 1,n&:j&=&Loop

            If j&=n&

                imprimir tab(j&*27);if(round(AI![i&,j&],8)<0,""," ");format$("%g",round(AI![i&,j&],8))

            Más

                imprimir tab(j&*27);if(round(AI![i&,j&],8)<0,""," ");format$("%g",round(AI![i&,j&],8)),

            endif

        endwhile

    endwhile

    'Calculate determinant
    S4000(n&,A![],D!)
    imprimir "\n\n Determinante = ";format$("%g",D!)
    'Check AP=A*A^-1=Unity Matrix
    imprimir "\n\n A*A^-1=AP=Einheismatrix:",
    ' AP=A*A^-1 matrix:'
    S3000(n&,A0![],AI![],AP![])

    whileloop 1,n&:i&=&Loop

        whileloop 1,n&:j&=&Loop

            If j&=n&

                imprimir tab(j&*27);if(round(AP![i&,j&],12)<0,""," ");format$("%g",round(AP![i&,j&],12))

            Más

                imprimir tab(j&*27);if(round(AP![i&,j&],12)<0,""," ");format$("%g",round(AP![i&,j&],12)),

            endif

        endwhile

    endwhile

ENDPROC'of main program

Proc S1000:parámetros n&,A![],V![]

    declarar S!,XA!,XB!,XG!,XL!
    imprimir

    whileloop 1,n&-1:k&=&Loop

        S! = 0

        whileloop k&,n&:i&=&Loop

            S!=S!+A![i&,K&]*A![i&,K&]

        endwhile

        XL!= ( (A![K&,K&]<0)-(A![K&,K&]>=0) ) *sqrt(S!)
        XA!= XL!*(XL!-A![K&,K&])
        V![K&]=A![K&,K&]-XL!

        whileloop k&+1,n&:i&=&Loop

            V![i&]=A![i&,K&]

        endwhile

        whileloop k&+1,2*n&:j&=&Loop

            XB!=0

            whileloop k&,n&:i&=&Loop

                XB!=XB!+V![i&]*A![i&,j&]

            endwhile

            XG!=XB!/XA!

            whileloop k&,n&:i&=&Loop

                A![i&,j&]=A![i&,j&]-XG!*V![i&]

            endwhile

        endwhile

        A![K&,K&]=XL!

        whileloop k&+1,n&:i&=&Loop

            A![i&,K&]=0

        endwhile

        imprimir "\n ";format$("0",K&);". Transformation:",

        whileloop 1,n&:i&=&Loop

            whileloop 1,n&:j&=&Loop

                If j&=n&

                    imprimir tab(j&*27);if(A![i&,j&]<0,""," ");format$("%g",A![i&,j&])

                Más

                    imprimir tab(j&*27);if(A![i&,j&]<0,""," ");format$("%g",A![i&,j&]),

                endif

            endwhile

        endwhile

    endwhile'k bucle

    volver

ENDPROC

Proc S2000 :parámetros n&,A![],B![],X![]

    declarar i&,j&,S!

    whileloop n&,1,-1:i&=&Loop

        j&=n&
        S!=0
        g20:
        caso i&=j&:goto "g40"
        S!=S!-A![i&,j&]*X![j&]
        j&=j&-1
        goto "g20"
        g40:
        X![i&] = (B![i&] + S!) / A![i&,i&]

    endwhile

    volver

ENDPROC

Proc S3000 :parámetros n&,A0![],AI![],AP![]

    declarar S!

    whileloop 1,n&:i&=&Loop

        whileloop 1,n&:j&=&Loop

            S!=0

            whileloop 1,n&:k&=&Loop

                S!=S!+A0![i&,k&]*AI![k&,j&]

            endwhile

            AP![i&,j&]=S!

        endwhile

    endwhile

    volver

ENDPROC

Proc S4000 :parámetros n&,A![],D!

    D! = 1

    whileloop 1,n&:i&=&Loop

        D!=D!*A![i&,i&]

    endwhile

    volver

ENDPROC