XProfan

Windowtitle "Newton-Kantorowitsch-Verfahren zur Lösung nichtlinearer Gleichungssysteme"
'{ (pD) private Demoware, 2013-03 von VB nach XProfan 11.2 by P. Specht
' Keine Gewähr! Original-Source: https://www.rhirte.de/vb/nlgleisys.htm
' ********************************************************************************************
' Kurzversion der Original-Erläuterung von Prof. Hirthe:
' Das Newtonverfahren zur Bestimmung einer Lösung (Nullstelle) kann für Systeme nichtlinearer
' Gleichungen verallgemeinert werden, indem man sich anstelle der einzelnen Funktionswerte
' Vektoren vorstellt. Statt der von Newton verwendeten 1. Ableitung der Funktion müssen dann
' aber vollständige partielle Ableitungen, zusammengefasst in der sog.  "Jacobi-Matrix" J[]
' bzw. der "Jacobi'schen Determinante" verwendet werden. Vorausgesetzt für n Variablen stehen
' n Bedingungsgleichungen zur Verfügung, wird die Lösung dann in allen n Dimensionen des
' Vektors gleichzeitig gesucht.
' Da man einen Vektor nicht durch eine Matrix dividieren kann, muß statt dessen mit ihrer
' Inversen multipliziert werden. Sei p[] der Vektor der x-Werte, F[] der Vektor der
' Funktionen, dann gilt ganz analog zur einfachen Newton-Formel die Iterationsformel:
' ********************************************************************************************
' *  p_next[1..n] = p[1..n] - J * (p[1..n])^-1 * F[p[1..n],1..n]
' ********************************************************************************************
' Wegen des Zeitaufwandes mag niemand die immer neuen Inversen berechnen,
' weshalb man diese in der Praxis nur ab und zu neu berechnet - das geht ja im einfachen
' Newtonverfahren mit der Ableitung auch (= temporärer Rückgriff auf regula falsi):
' Wenn man der Differenz von p und p_next den Namen "Verbesserungsvektor u[]" gibt und in
' der Formel alles von links mit J[] multipliziert, so ergibt J*J^-1 eine Einheitsmatrix, die
' man weglassen kann. Man erhält: Ein LINEARES Gleichungssystem in 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]>] ]
' ********************************************************************************************
' Man kann nichtlineare Probleme also linearisieren: u[] stabilisiert sich in Nähe der Lösung!
' ********************************************************************************************
' Anmerkung: "Ableitung" ist als numerische Ableitung programmiert, weil man damit unabhängig
' vom konkreten Funktionensatz ist, dessen Lösung gefunden werden soll. Anders "Gleichungen",
' hier stehen die veränderlichen Gleichungen, die man leider nur unter großen Schwierigkeiten
' zur Laufzeit des Programmes eintragen könnte.
'
' Für die Lineare Lösung wird dann ein übliches Gauss-Jordan- oder Gauss-Seidel-Verfahren benutzt.
'}  u[] ist stabil, wenn der Unterschied z.B. kleiner eps = 1E-9 beträgt.
'{ Deklarationen, Initialisierung
' für die Iteration x_next[]=J*x[]
Font 2:set("decimals",15):randomize:cls rnd(8^8):window 0,0-%maxx,480
Declare n&,i&,j&,x$[],Runde&,eps!,summe!
Init:
x$[]=explode("0.1,0.1,0.1",",")
n&=sizeof(x$[]):declare B![n&,n&+1],x![n&],u![n&]

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

    eps!=10^-9
    '}
    '{ Hauptteil mit Linearisierung und numerischer Berechnung der Jacobi-Matrix
    Runde&=0

    Repeat

        whileloop n&:i&=&Loop:whileloop n&:j&=&Loop

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

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

            Until (Summe!<eps!) OR (Runde&>=2000)' Loesung_ausgeben

            Case Runde&>=2000:Print "Auch nach 2000 Iterationen kein stabiles Ergebnis!"

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

                Waitinput
                '}
                End

                Proc Gleichung :Parameters num&,x![]' Die n& Gleichungen hier eingeben:

                    Declare g!' ACHTUNG: Sqr = Quadrat (Wurzel wäre sqrt()!)

                    Select num&' Probe: Die 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
                        Print " Error: Eine Gleichung ";num&,"gibt es nicht!":waitinput :waitinput :end

                    EndSelect

                    ' print num&,x![1],x![2],x![3],g! ' DEBUG Line
                    Return g!

                Endproc

                Proc Ableitung :Parameters n&,Gleinum&,Varnum&,xv![]

                    Declare dx!,y1!,y2!,dxv![10],i&

                    Whileloop n&:i&=&Loop

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

                    Endwhile

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

                Endproc

                Proc GaussPivotElimination :Parameters n%,a![],x![]'Eines der hier möglichen Gauss-...-Verfahren

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

                    WhileLoop n%:i%=&Loop:merk%[i%]=i%:EndWhile

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

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

                            EndWhile

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

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

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

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

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

                                        EndIf

                                    EndWhile

                                EndWhile

                                If jmax%<>k%

                                    WhileLoop k%,n%+1:j%=&Loop: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

                                    If kmax% <> k%

                                        WhileLoop n%:i%=&Loop: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%=&Loop:s!=A![i%,k%]/A![k%,k%]:A![i%,k%]=0

                                            WhileLoop k%+1,n%+1:j%=&Loop: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%=&Loop:s! = A![i%,n%+1]

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

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

                                            EndWhile

                                            return x![]

                                        endproc