XProfan

Windowtitle "Newton-Kantorowitsch-take action to Solution nichtlinearer Gleichungssysteme"
'{ (pD) private Demoware, 2013-03 of VB to XProfan 11.2 by P. woodpecker
' No Gewähr! Original-Source: https://www.rhirte.de/vb/nlgleisys.htm
' ********************************************************************************************
' Kurzversion the Original-explanation of Prof. Hirthe:
' the Newtonverfahren to Bestimmung of/ one Solution (Nullstelle) can for systems nichtlinearer
' Gleichungen generalizing go, because one itself in lieu of the individual Funktionswerte
' Vektoren vorstellt. but not the of Newton verwendeten 1. Ableitung the function must then
' but complete partial Ableitungen, abstracted in the undertow.  "Jacobi-Matrix" J[]
' or. the "Jacobi'schen Determinante" uses go. presupposed for n variables stand
' n Bedingungsgleichungen available, becomes The Solution then in all n Dimensionen the
' Vektors at the same time sought.
' there one a vector not by a Matrix divide can, must instead with yours
' Inversen multipliziert go. be p[] the vector the x-values, F[] the vector the
' functions, then counts integrally analog to einfachen Newton-Formel The Iterationsformel:
' ********************************************************************************************
' *  p_next[1..n] = p[1..n] - J * (p[1..n])^-1 * F[p[1..n],1..n]
' ********************************************************************************************
' because of the Zeitaufwandes likes nobody the ever new Inversen to charge,
' and so one these in the practice only now and than recalculated - the goes Yes in the einfachen
' Newtonverfahren with the Ableitung too (= temporärer Rückgriff on regula falsi):
' If one the difference of p and p_next whom names "Verbesserungsvektor u[]" gives and
' the Formel everything of left with J[] multipliziert, so yields J*J^-1 a Einheitsmatrix, The
' one omit can. one sustain: One LINEARES Gleichungssystem in vektorieller spelling!
'   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]>] ]
' ********************************************************************************************
' one can nichtlineare Problems means linearisieren: u[] stabilisiert itself in sew the Solution!
' ********************************************************************************************
' Note: "Ableitung" is as numerische Ableitung program, because one so independent
' of concrete Funktionensatz is, which Solution found go should. differently "Gleichungen",
' here stand The veränderlichen Gleichungen, everybody can unfortunately only under large difficulty
' to Laufzeit the Program present could.
'
' for the Lineare Solution becomes then one übliches Gauss-Jordan- or Gauss-Seidel-take action using.
'}  u[] is steady, if the Difference z.B. small eps = 1E-9 totals.
'{ Deklarationen, Initialisierung
' for 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!,amount!
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 with Linearisierung and numerischer Berechnung the 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' one the linearen Gauß-...-take action
            GaussPivotElimination n&,B![],u![]
            amount!=0

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

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

            Case Runde&>=2000:Print "Auch to 2000 Iterationen no stabiles Result!"

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

                Waitinput
                '}
                End

                Proc Gleichung :Parameters num&,x![]' The n& Gleichungen here prompt:

                    Declare g!' ACHTUNG: Sqr = Quadrat (root would sqrt()!)

                    Select num&' Probe: The Musterlösung must x[1]=1, x[2]=-2, x[3]=4 yield!

                        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: an Gleichung ";num&,"gibt not!":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![]'one the here possible Gauss-...-take action

                    ' a![n%,(n%+1)] ' LGS-Lines*(Split+1) wg.Lösungsvektor "Rechte Page "
                    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