XProfan

Français

p.specht

qui Householder-Algorithmus
=====================
chez kritischen Datenkonstellationen versagen übliche Lösungsverfahren pour Gleichungssysteme. cela Gauss'sche Eliminationsverfahren etwa produziert chez grand Werte-Spannweiten dedans qui Matrizenspalten erhebliche faute, qui trop absoluter Unbrauchbarkeit qui Ergebnisse mener peut, sans dass qui Anwender es juste merken serait. Pour solche Fälle hatte qui États-Unis-Amerikanische Mathematiker Alston Householder 1958 qui concept, qui pas-Diagonalwerte par envoyé gewählte Spiegelungen trop zéro trop faire. Besonders effizient gelingt cela chez Diagonalsymmetrischen Matrizen, comment vous dans qui Physik souvent vorkommen. (Ansonsten peux on chacun beliebige Matrix par Multiplikation avec ihrer Transponierten trop einer Diagonalsymmetrischen Matrix faire et cela Resultat ensuite LU-décomposer).

dans cette Demo pour private Bildungszwecke wurde Fortran's Arrayindex-la base '1' beibehalten! dans XProfan une ziemliche Stolperfalle!

Kompilieren Marque Séparation

Titre de la fenêtre "Inversion einer rellen, diagonalsymetrischen Matrix mittels Householder-Algorithmus"
'Solution per 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 pour XProfan11.2a 2014-10 by P.Specht, Wien;  OHNE JEDWEDE GEWÄHR!
'Verfahrensdetails siehe https://de.wikipedia.org/wiki/Householdertransformation
'***************************************************************
'{  Inversion of a réel square matrix by Householder's method  *
'
'* SAMPLE RUN [im Original seulement single precision, Anm. P.Specht]*
'* 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:
'* Contribution 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 Version By J-P Moreau, Paris.   *
'*                            (www.jpmoreau.fr)                *
'*       XProfan11.2a Migration by P.Specht, Vienna            *
'}
Fenêtre Style 24
Fenêtre 0,0-%maxx,%maxy-40
Font 2
set("decimals",17)
Déclarer A![25,50]
Imprimer "\n\n               DER HOUSHOLDER-ALGORITHMUS\n               =========================="
Householder(0,A![])' chez Angabe n = 0 Dimensionen wird cela Testbeispiel berechnet
Waitinput
FIN

Proc Householder :parameters n&,A![]

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

    si n&=0'Read matrix to être 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

    Imprimer "\n ist pour diagonalsymetrische Ausgangsmatrizen vorgesehen,   BEISPIEL:\n"

    whileloop 1,n&:i&=&Boucle

        whileloop 1,n&:j&=&Boucle

            Si j& = n&

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

            D'autre

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

            endif

        endwhile

    endwhile

    whileloop 1,n&:i&=&Boucle

        whileloop 1,n&:j&=&Boucle

            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&=&Boucle

        whileloop 1,n&:i&=&Boucle

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

        endwhile

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

        whileloop 1,n&:i&=&Boucle

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

        endwhile

    endwhile'k loop

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

    whileloop 1,n&:i&=&Boucle

        whileloop 1,n&:j&=&Boucle

            Si j&=n&

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

            D'autre

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

            endif

        endwhile

    endwhile

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

    whileloop 1,n&:i&=&Boucle

        whileloop 1,n&:j&=&Boucle

            Si j&=n&

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

            D'autre

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

            endif

        endwhile

    endwhile

ENDPROC'of main program

Proc S1000:parameters n&,A![],V![]

    declare S!,XA!,XB!,XG!,XL!
    imprimer

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

        S! = 0

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

            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&=&Boucle

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

        endwhile

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

            XB!=0

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

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

            endwhile

            XG!=XB!/XA!

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

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

            endwhile

        endwhile

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

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

            A![i&,K&]=0

        endwhile

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

        whileloop 1,n&:i&=&Boucle

            whileloop 1,n&:j&=&Boucle

                Si j&=n&

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

                D'autre

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

                endif

            endwhile

        endwhile

    endwhile'k loop

    return

ENDPROC

Proc S2000 :parameters n&,A![],B![],X![]

    declare i&,j&,S!

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

        j&=n&
        S!=0
        g20:
        cas 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

    return

ENDPROC

Proc S3000 :parameters n&,A0![],AI![],AP![]

    declare S!

    whileloop 1,n&:i&=&Boucle

        whileloop 1,n&:j&=&Boucle

            S!=0

            whileloop 1,n&:k&=&Boucle

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

            endwhile

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

        endwhile

    endwhile

    return

ENDPROC

Proc S4000 :parameters n&,A![],D!

    D! = 1

    whileloop 1,n&:i&=&Boucle

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

    endwhile

    return

ENDPROC