XProfan

Windowtitle "Newton-Kantorowitsch-Verfahren a Solución nichtlinearer Gleichungssysteme"
'{ (pD) private Demoware, 2013-03 de VB después de XProfan 11.2 by P. Pájaro carpintero
' Keine Gewähr! Original-Source: https://www.rhirte.de/vb/nlgleisys.htm
' ********************************************************************************************
' Kurzversion el Original-Erläuterung de Prof. Hirthe:
' Das Newtonverfahren a Bestimmung uno Solución (Nullstelle) kann para Systeme nichtlinearer
' Gleichungen verallgemeinert voluntad, indem uno se anstelle el individual Funktionswerte
' Vektoren vorstellt. Statt el de Newton verwendeten 1. Ableitung el Función necesario entonces
' aber vollständige partielle Ableitungen, zusammengefasst en el sog.  "Jacobi-Matrix" J[]
' o. el "Jacobi'schen Determinante" verwendet voluntad. Vorausgesetzt para n Variables posición
' n Bedingungsgleichungen disponible, se el Solución entonces en allen n Dimensionen des
' Vektors gleichzeitig gesucht.
' Como uno una Vektor no por una Matrix dividieren kann, muß en lugar de dessen con ihrer
' Inversen multipliziert voluntad. Sei p[] el Vektor el x-Werte, F[] el Vektor el
' Características, entonces gilt bastante analog a einfachen Newton-Formel el Iterationsformel:
' ********************************************************************************************
' *  p_next[1..n] = p[1..n] - J * (p[1..n])^-1 * F[p[1..n],1..n]
' ********************************************************************************************
' Wegen des Zeitaufwandes mag niemand el siempre neuen Inversen berechnen,
' por qué uno esta en el Praxis sólo de y a neu berechnet - el va sí en el einfachen
' Newtonverfahren con el Ableitung auch (= temporärer Rückgriff en regula falsi):
' Wenn uno el Differenz de p y p_next el Namen "Verbesserungsvektor u[]" son y
' el Formel alles de links con J[] multipliziert, así ergibt J*J^-1 una Einheitsmatrix, el
' uno weglassen kann. Man erhält: Ein LINEARES Gleichungssystem en 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]>] ]
' ********************************************************************************************
' Se puede nichtlineare Problemas also linearisieren: u[] stabilisiert se en Nähe el Solución!
' ********************************************************************************************
' Anmerkung: "Ableitung" es como numerische Ableitung programmiert, porque uno así unabhängig
' vom konkreten Funktionensatz es, dessen Solución gefunden voluntad se. Anders "Gleichungen",
' hier posición el veränderlichen Gleichungen, el uno desafortunadamente sólo bajo großen Schwierigkeiten
' a Laufzeit des Programmes eintragen podría.
'
' Für el Lineare Solución se entonces una übliches Gauss-Jordan- oder Gauss-Seidel-Verfahren benutzt.
'}  u[] es stabil, si la Diferencia z.B. kleiner eps = 1E-9 beträgt.
'{ Deklarationen, Initialisierung
' para el Iteration x_next[]=J*x[]
Font 2:set("decimals",15):randomize:cls rnd(8^8):window 0,0-%maxx,480
Declarar n&,i&,j&,x$[],Runde&,eps!,summe!
Init:
x$[]=explode("0.1,0.1,0.1",",")
n&=sizeof(x$[]):declarar B![n&,n&+1],x![n&],u![n&]

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

    eps!=10^-9
    '}
    '{ Hauptteil con Linearisierung y numerischer Berechnung el 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 el 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!) O (Runde&>=2000)' Loesung_ausgeben

            Case Runde&>=2000:Imprimir "Auch después de 2000 Iterationen kein stabiles Ergebnis!"

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

                Waitinput
                '}
                End

                Proc Gleichung :Parámetros num&,x![]' El n& Gleichungen hier eingeben:

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

                    Select num&' Probe: El 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
                        Imprimir " Error: Un Gleichung ";num&,"gibt no!":waitinput :waitinput :end

                    EndSelect

                    ' imprimir num&,x![1],x![2],x![3],g! ' DEBUG Line
                    Volver g!

                ENDPROC

                Proc Ableitung :Parámetros n&,Gleinum&,Varnum&,xv![]

                    Declarar 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![])
                    Volver (y2!-y1!)/dx!

                ENDPROC

                Proc GaussPivotElimination :Parámetros n%,a![],x![]'Eines el hier möglichen Gauss-...-Verfahren

                    ' a![n%,(n%+1)] ' LGS-Zeilen*(Spalten+1) ca.Lösungsvektor "Rechte Seite"
                    declarar 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

                                            volver x![]

                                        ENDPROC