XProfan

Windowtitle "Newton-Kantorowitsch-procéder zur Solution nichtlinearer Gleichungssysteme"
'{ (pD) private Demoware, 2013-03 de VB pour XProfan 11.2 by P. Specht
' aucun Gewähr! Original-Source: https://www.rhirte.de/vb/nlgleisys.htm
' ********************************************************************************************
' Kurzversion qui Original-Erläuterung de Prof. Hirthe:
' cela Newtonverfahren zur Bestimmung einer Solution (Nullstelle) peux pour Systeme nichtlinearer
' Gleichungen verallgemeinert volonté, indem on sich anstelle qui individuel Funktionswerte
' Vektoren vorstellt. Statt qui de Newton verwendeten 1. Ableitung qui Funktion doit ensuite
' mais vollständige partielle Ableitungen, zusammengefasst dans qui sog.  "Jacobi-Matrix" J[]
' bzw. qui "Jacobi'schen Determinante" verwendet volonté. Vorausgesetzt pour n Variablen stehen
' n Bedingungsgleichungen zur Disposition, wird qui Solution ensuite dans allen n Dimensionen des
' Vektors gleichzeitig gesucht.
' là on une Vektor pas par une Matrix dividieren peux, doit statt dessen avec ihrer
' Inversen multipliziert volonté. Sei p[] qui Vektor qui x-Werte, F[] qui Vektor qui
' Funktionen, ensuite gilt entier analog zur einfachen Newton-Formel qui Iterationsformel:
' ********************************************************************************************
' *  p_next[1..n] = p[1..n] - J * (p[1..n])^-1 * F[p[1..n],1..n]
' ********************************************************************************************
' à cause de des Zeitaufwandes mag niemand qui toujours neuen Inversen berechnen,
' weshalb on cet dans qui Praxis seulement ab et trop récente berechnet - cela allez oui im einfachen
' Newtonverfahren avec qui Ableitung aussi (= temporärer Rückgriff sur regula falsi):
' si on qui Differenz de p et p_next den Namen "Verbesserungsvektor u[]" gibt et dans
' qui Formel alles de à gauche avec J[] multipliziert, so ergibt J*J^-1 une Einheitsmatrix, qui
' on omettre peux. on erhält: un LINEARES Gleichungssystem dans vektorieller Schreibweise!
'   J[ <delta(F[p[i]])/delta(x[j=1..n])>, x[j=1..n] ] * u[i=1..n] = -1 * F[ p[<x_[i=1..n]>] ]
' ********************************************************************************************
' il peut nichtlineare Probleme alors linearisieren: u[] stabilisiert sich dans Nähe qui Solution!
' ********************************************************************************************
' Anmerkung: "Ableitung" ist comme numerische Ableitung programmiert, weil on avec cela indépendant
' vom konkreten Funktionensatz ist, dessen Solution trouvé volonté soll. Anders "Gleichungen",
' ici stehen qui veränderlichen Gleichungen, qui on malheureusement seulement sous grand Schwierigkeiten
' zur Laufzeit des Programmes eintragen pourrait.
'
' Pour qui Lineare Solution wird ensuite un übliches Gauss-Jordan- ou bien Gauss-Seidel-procéder benutzt.
'}  u[] ist stabil, si qui Unterschied z.B. kleiner eps = 1E-9 beträgt.
'{ Deklarationen, initialisation
' pour qui Iteration x_next[]=J*x[]
Font 2:set("decimals",15):randomize:cls rnd(8^8):window 0,0-%maxx,480
Déclarer n&,i&,j&,x$[],Runde&,eps!,somme!
Init:
x$[]=explode("0.1,0.1,0.1",»)
n&=sizeof(x$[]):declare B![n&,n&+1],x![n&],u![n&]

whileloop n&:x![&Boucle]=val(x$[&Boucle-1]):endwhile :clear x$[]

    eps!=10^-9
    '}
    '{ Hauptteil avec Linearisierung et numerischer Berechnung qui Jacobi-Matrix
    Runde&=0

    Repeat

        whileloop n&:i&=&Boucle:whileloop n&:j&=&Boucle

            B![i&,j&]=Ableitung(n&,i&,j&,x![]):endwhile
            B![i&,n&+1]= -1*Gleichung(i&,x![]):endwhile' Eines qui linearen Gauß-...-procéder
            GaussPivotElimination n&,B![],u![]
            somme!=0

            whileloop n&:i&=&Boucle:x![i&]=u![i&]+x![i&]:somme!=somme!+Abs(u![i&]):endwhile :inc Runde&

            Until (somme!<eps!) OU (Runde&>=2000)' Loesung_ausgeben

            Cas Runde&>=2000:Imprimer "Auch pour 2000 Iterationen ne...aucune stabiles Ergebnis!"

            Whileloop n&:i&=&Boucle:imprimer "x"+str$(i&);" = ";format$("%g",x![i&]):endwhile

                Waitinput
                '}
                Fin

                Proc Gleichung :Paramètres num&,x![]' qui n& Gleichungen ici eingeben:

                    Déclarer g!' ACHTUNG: Sqr = Quadrat (Wurzel wäre sqrt()!)

                    Select num&' Probe: qui Musterlösung muss x[1]=1, x[2]=-2, x[3]=4 ergeben!

                        Caseof 1

                        g!=3*x![1]+4*sqr(x![2])-6*x![3]+5

                        Caseof 2

                        g!=sqr(x![1])-3*x![2]+5*x![3]-27

                        Caseof 3

                        g!= -5*x![1]+x![2]+sqr(x![3])- 9
                        Otherwise
                        Imprimer " Error: une Gleichung ";num&,"gibt es pas!":waitinput :waitinput :end

                    EndSelect

                    ' imprimer num&,x![1],x![2],x![3],g! ' DEBUG Line
                    Retour g!

                ENDPROC

                Proc Ableitung :Paramètres n&,Gleinum&,Varnum&,xv![]

                    Déclarer dx!,y1!,y2!,dxv![10],i&

                    Whileloop n&:i&=&Boucle

                        dxv![i&] = xv![i&]

                    Endwhile

                    y1! = Gleichung(Gleinum&,xv![])
                    dx! = eps! * xv![Varnum&]
                    Cas dx!=0:dx!=eps!
                    dxv![Varnum&] = dxv![Varnum&] + dx!
                    y2!=Gleichung(Gleinum&,dxv![])
                    Retour (y2!-y1!)/dx!

                ENDPROC

                Proc GaussPivotElimination :Paramètres n%,a![],x![]'Eines qui ici möglichen Gauss-...-procéder

                    ' a![n%,(n%+1)] ' LGS-Zeilen*(Spalten+1) GT.Lösungsvektor «Les droits Seite"
                    declare i%,j%,k%,jmax%,kmax!,kmax%,merk%[n%],s!,max!,skal![n%]

                    WhileLoop n%:i%=&Boucle:merk%[i%]=i%:Endwhile

                        WhileLoop n%:i%=&Boucle:s!=0

                            WhileLoop n%:j%=&Boucle:s!=s!+Abs(A![i%,j%]):Endwhile :skal![i%]=1/s!

                            Endwhile

                            WhileLoop n%-1:k%=&Boucle:max!=skal![k%]*Abs(A![k%,k%]):kmax%=k%:jmax%=k%

                                WhileLoop k%,n%:j%=&Boucle

                                    WhileLoop k%,n%:i%=&Boucle

                                        Si (skal![j%]*Abs(A![j%,i%]))>max!

                                            jmax%=j%:kmax%=i%:max!=skal![j%]*Abs(A![j%,i%])

                                        EndIf

                                    Endwhile

                                Endwhile

                                Si jmax%<>k%

                                    WhileLoop k%,n%+1:j%=&Boucle:s!=A![k%,j%]:A![k%,j%]=A![jmax%,j%]:A![jmax%,j%]=s!

                                        Endwhile :s!=skal![k%]:skal![k%]=skal![jmax%]:skal![jmax%] = s!

                                    EndIf

                                    Si kmax% <> k%

                                        WhileLoop n%:i%=&Boucle:s!=A![i%,k%]:A![i%,k%]=A![i%,kmax%]:A![i%,kmax%]=s!

                                            Endwhile :j% = merk%[k%]:merk%[k%]=merk%[kmax%]:merk%[kmax%]=j%

                                        EndIf

                                        WhileLoop k%+1,n%:i%=&Boucle:s!=A![i%,k%]/A![k%,k%]:A![i%,k%]=0

                                            WhileLoop k%+1,n%+1:j%=&Boucle:A![i%,j%]=A![i%,j%]-s!*A![k%,j%]:Endwhile

                                            Endwhile

                                        Endwhile

                                        x![merk%[n%]]=A![n%,n%+1]/A![n%,n%]

                                        WhileLoop n%-1,1,-1:i%=&Boucle:s! = A![i%,n%+1]

                                            WhileLoop i%+1,n%:j%=&Boucle:s!=s!-A![i%,j%]*x![merk%[j%]]:Endwhile

                                                x![merk%[i%]]=s!/A![i%,i%]

                                            Endwhile

                                            return x![]

                                        endproc