XProfan

Deutsch

p.specht

Bei der rechnerischen Lösung einer u.U. sehr komplizierten Gleichung 'F(x)=Vorgabewert', also der Suche all jener x-Werte, bei denen diese Gleichung erfüllt bzw. 'wahr' ist, ordnet man üblicherweise die Gleichung um zu f(x) minus Vorgabewert = 0, und sucht nach den Nullstellen des y-Funktionswertes. Das BISECTION(Unterteilungs)-Verfahren war dabei einer der allerersten erfolgreich angewendeten Computeralgorithmen, und hat sich bis heute als sehr robuste Methode bewährt, soll heißen: Das Verfahren liefert - abhängig vom angebenen Suchbereich - mit Sicherheit genau einen exakten Nulldurchgang der Funktion als Lösung, falls eine solche Lösung existiert.

Weiterer Vorteil: Die gefundenen Lösungen gelten als "technisch stabil". Diese Stabilität wird aus Sicht der Mathematik-Puristen allerdings duch einen Nachteil erkauft: Es werden nur Lösungen gefunden, die die x-Achse tatsächlich schneiden, also in der Nähe der Lösungspunkte sowohl positive als auch negative y-Achsenwerte besitzen. Mit anderen Worten: Der Algorithmus kann keine Lösungen ermitteln, die durch bloße Berührung der x-Achse nur von oben oder nur von unten gegeben wären. In solchen Fällen wäre Verfahren wie "Regula falsi" oder "Newton-Raphson" geeigneter, beide haben allerdings auch wieder ihre Vor- und Nachteile.

Hinweis: Die Urheberrechte liegen bei ACM; hier wurde nur die Portierbarkeit nach XProfan demonstriert. Ohne jede Gewähr daher:

Kompilieren Markieren Separieren

' ACM TOMS Algorithm 4 'BISECTION'
'{Depending on 4 start parameters, BISEC finds one x where Function y=F(x)=0
' Original Fortran Source by S. Gorn from https://www.netlib.org/toms/4.gz
' (C) Feb.1960 by ACM-TOMS Association of Computing Machinery
' ----------------------------------------------------------------------------
' (D) Demo-Migration to XProfan-11.2a in 2014-11 by P.Specht, Vienna (Austria)
' Es bestehen Urheberrechte Dritter! No Warranty Whatsoever! Ohne jede Gewähr!
' ----------------------------------------------------------------------------
' Description (F99): Date: Sun, 12 Feb 95 18:54:30 +0000
' Here is a transcription of algorithm #4 from Collected Algorithms
' from ACM. I had to modify some characters to make it fit into ASCII,
' specially the multiplicative operator has become * and greek letters
' are substituted by their equivalent names.
'         Jose R. Valverde, European Bioinformatics Institute, txomsy@ebi.ac.uk
' ------------------------------------------------------------------------------
' Original Algorithm 4: BISECTION ROUTINE
' by S. Gorn, Univeristy of Pennsylvania Computer Center, Philadelphia, Pa.
' ------------------------------------------------------------------------------
' Comment: This procedure evaluates a function at the end-points
' 	of a real interval, switching to an error exit (fools
' 	exit) FLSXT if there is no change of sign. Otherwise
' 	it finds a root by iterated bisection and evaluation
' 	at the midpoint, halting if either the value of the
' 	function is less than the free variable epsilon, or two
' 	successive approximations of the root differ by less
' 	than epsilon1. Epsilon should be chosen of the order of
' 	error in evaluating the function (otherwise time would be
' 	wasted), and epsilon1 of the order of desired accuracy.
' 	Epsilon1 must not be less than two units in the last place
' 	carried by the machine or else indefinite cycling will
' 	occur due to roundoff on bisection. Although this
' 	method is of 0 order, and therefore among the slow-
' 	est, it is applicable to any continuous function. The
' 	fact that no differentiability conditions have to be
' 	checked makes it, therefore, an 'old work-horse'
' 	among routines for finding real roots which have
' 	already been isolated. The free varaibles y1 and y2
' 	are (presumably) the end-points of an interval within
' 	which there is an odd number of roots of the real function F.
' 	Alpha is the temporary exit fot the evaluation of F.;
'---------------------------------------------------------------
'---------------------------------------------------------------
' CERTIFICATION OF ALGORITHM 4
' BISECTION ROUTINE (S. Gorn, _Comm._ACM_, March 1960), Program by
' Patty Jane Rader,* Argonne National Laboratory, Argonne, Illinois
'
' Bisec was coded for the Royal Precision LGP-30 computer, using inter-
' pretive floating point system (24.2) with 28 bits of significance.
'
'  The following minor correction was found necessary.
'	alpha: go to gamma[1] should be go to gamma[i]
' After this correction was made, the program ran smoothly for
'  F(x) = cos x, using the following parameters:
'---------------------------------------------------------------
'	y1	 y2	 Epsilon	Epsilon1      Results
'---------------------------------------------------------------
' 0     1	  .001		.001		      FLSXT
' 0	    2	  .001		.001		      1.5703
'	1.5	  2	  .001		.001		      1.5703
'	1.55	2	  .1		  .1		        1.5500
'	1.5	  2	  .001		.1		        1.5625
'---------------------------------------------------------------
' These combinations test all loops of the program.
' *) Work supported by the U. S. Atomic Energy Commission.
'}
'---------------------------------------------------------------
'{TESTPROGRAMM
CLS
Font 2
Set("decimals",15)
print "\n Berechnete Testwerte bitte mit den im Kommentar"
print " angegebenen Tabellenresultaten vergleichen: \n"
print "   ";format$("%g",BiSec(0,1,0.001,0.001))
print "   ";format$("%g",BiSec(0,2,0.001,0.001))
print "   ";format$("%g",BiSec(1.5,2,0.001,0.001))
print "   ";format$("%g",BiSec(1.55,2,0.1,0.1))
print "   ";format$("%g",BiSec(1.5,2,0.001,0.1))
Waitinput
END
'}-------------------------------------------------------------------------
'Function F(x) HIER EINPORGRAMMIEREN!

proc F :parameters x!

    return cos(x!)

endproc

'-------------------------------------------------------------------------
'Bisec Routine:

proc BISEC :parameters y1!,y2!,epsilon!,epsilon1!

    declare x!,f!,f1!,i&,j&,k&,FLSXT!
    FLSXT!=-999999999'bedeutet ERROR, keine Nullstelle in diesem Bereich
    Bisec:
    i&=1:j&=1:k&=1:x!=y2!
    Alpha:
    f!=F(x!) : case abs(f!)<=epsilon!:return x!

    if i&=1:goto "First_val" : elseif i&=2:goto "Succ_val" :endif

        First_val:
        i&=2 : f1!=f! : x!=y1!
        goto "Alpha"
        Succ_val:

        if (f!*f1!)>=0

            if j&=1: Print "Keine Nullstelle in diesem Bereich! Error ";:return FLSXT!

            elseif j&=2: goto "Reg_delta"

            endif

        endif

        if k&=1:goto "Sec_val" : elseif k&=2:goto "Reg_etha" :endif

            Sec_val:
            j&=2:k&=2
            Midpoint:
            x!=(y1!+y2!)/2
            goto "Alpha"
            Reg_delta:
            y2!=x!
            Precision:
            case (abs(y1!-y2!)>=epsilon1!):goto "Midpoint"
            return x!
            Reg_etha:
            y1!=x!
            goto "Precision"

        EndProc

XProfan 11
Computer: Gerät, daß es in Mikrosekunden erlaubt, 50.000 Fehler zu machen, zB 'daß' statt 'das'...

16.05.2021 ▲

Vollständiges Zitat

Absatzweise zitieren Zeilenweise zitieren

p.specht

Michael Wodrich hat eine verbesserte Variante vorgestellt:

Sein Kommentar: "... und ein zweiter BiSec(t) ohne Goto aber mit Iterationen-Begrenzer (damit es nicht endlos läuft). Beides zusammen gemixt und eine Iteration vorher ausgebremst..."

Link: [...]

Kompilieren Markieren Separieren

' ACM TOMS Algorithm 4 'BISECTION'
'{Depending on 4 start parameters, BISEC finds one x where Function y=F(x)=0
' Original Fortran Source by S. Gorn from https://www.netlib.org/toms/4.gz
' (C) Feb.1960 by ACM-TOMS Association of Computing Machinery
' ----------------------------------------------------------------------------
' (D) Demo-Migration to XProfan-11.2a in 2014-11 by P.Specht, Vienna (Austria)
' Es bestehen Urheberrechte Dritter! No Warranty Whatsoever! Ohne jede Gewähr!
' ----------------------------------------------------------------------------
' Description (F99): Date: Sun, 12 Feb 95 18:54:30 +0000
' Here is a transcription of algorithm #4 from Collected Algorithms
' from ACM. I had to modify some characters to make it fit into ASCII,
' specially the multiplicative operator has become * and greek letters
' are substituted by their equivalent names.
'         Jose R. Valverde, European Bioinformatics Institute, txomsy@ebi.ac.uk
' ------------------------------------------------------------------------------
' Original Algorithm 4: BISECTION ROUTINE
' by S. Gorn, Univeristy of Pennsylvania Computer Center, Philadelphia, Pa.
' ------------------------------------------------------------------------------
' Comment: This procedure evaluates a function at the end-points
' 	of a real interval, switching to an error exit (fools
' 	exit) FLSXT if there is no change of sign. Otherwise
' 	it finds a root by iterated bisection and evaluation
' 	at the midpoint, halting if either the value of the
' 	function is less than the free variable epsilon, or two
' 	successive approximations of the root differ by less
' 	than epsilon1. Epsilon should be chosen of the order of
' 	error in evaluating the function (otherwise time would be
' 	wasted), and epsilon1 of the order of desired accuracy.
' 	Epsilon1 must not be less than two units in the last place
' 	carried by the machine or else indefinite cycling will
' 	occur due to roundoff on bisection. Although this
' 	method is of 0 order, and therefore among the slow-
' 	est, it is applicable to any continuous function. The
' 	fact that no differentiability conditions have to be
' 	checked makes it, therefore, an 'old work-horse'
' 	among routines for finding real roots which have
' 	already been isolated. The free varaibles y1 and y2
' 	are (presumably) the end-points of an interval within
' 	which there is an odd number of roots of the real function F.
' 	Alpha is the temporary exit fot the evaluation of F.;
'---------------------------------------------------------------
'---------------------------------------------------------------
' CERTIFICATION OF ALGORITHM 4
' BISECTION ROUTINE (S. Gorn, _Comm._ACM_, March 1960), Program by
' Patty Jane Rader,* Argonne National Laboratory, Argonne, Illinois
'
' Bisec was coded for the Royal Precision LGP-30 computer, using inter-
' pretive floating point system (24.2) with 28 bits of significance.
'
'  The following minor correction was found necessary.
'	alpha: go to gamma[1] should be go to gamma[i]
' After this correction was made, the program ran smoothly for
'  F(x) = cos x, using the following parameters:
'---------------------------------------------------------------
'	y1	 y2	 Epsilon	Epsilon1      Results
'---------------------------------------------------------------
' 0     1	  .001		.001		      FLSXT
' 0	    2	  .001		.001		      1.5703
'	1.5	  2	  .001		.001		      1.5703
'	1.55	2	  .1		  .1		        1.5500
'	1.5	  2	  .001		.1		        1.5625
'---------------------------------------------------------------
' These combinations test all loops of the program.
' *) Work supported by the U. S. Atomic Energy Commission.
'}
'---------------------------------------------------------------
'{TESTPROGRAMM
CLS
Font 2
Set("decimals",15)
print "---------------------------------------------------------------"
print Tab(2);"y1";  Tab(8);"y2";Tab(12);"Epsilon";Tab(21);"Epsilon1";Tab(31);"Results"
print "---------------------------------------------------------------"
print Tab(2);"0";   Tab(8);"1"; Tab(12);".001";   Tab(21);".001";    Tab(31);"FLSXT"
print Tab(2);"0";   Tab(8);"2"; Tab(12);".001";   Tab(21);".001";    Tab(31);"1.5703"
print Tab(2);"1.5"; Tab(8);"2"; Tab(12);".001";   Tab(21);".001";    Tab(31);"1.5703"
print Tab(2);"1.55";Tab(8);"2"; Tab(12);".1";     Tab(21);".1";      Tab(31);"1.5500"
print Tab(2);"1.5"; Tab(8);"2"; Tab(12);".001";   Tab(21);".1";      Tab(31);"1.5625"
print "---------------------------------------------------------------"
print "\n Berechnete Testwerte bitte mit den obigen"
print " Tabellenresultaten vergleichen: \n"
print format$("   BiSec(0,    1, .001, .001) = %g",BiSec(0,1,0.001,0.001))
print format$("   BiSec(0,    2, .001, .001) = %g",BiSec(0,2,0.001,0.001))
print format$("   BiSec(1.5,  2, .001, .001) = %g",BiSec(1.5,2,0.001,0.001))
print format$("   BiSec(1.55, 2, .1,   .1  ) = %g",BiSec(1.55,2,0.1,0.1))
print format$("   BiSec(1.5,  2, .001, .1  ) = %g",BiSec(1.5,2,0.001,0.1))
print "---------------------------------------------------------------"
print format$("   BiSect(0,    1, .001, .001) = %g",BiSect(0,1,0.001,0.001))
print format$("   BiSect(0,    2, .001, .001) = %g",BiSect(0,2,0.001,0.001,9))
print format$("   BiSect(1.5,  2, .001, .001) = %g",BiSect(1.5,2,0.001,0.001))
print format$("   BiSect(1.55, 2, .1,   .1  ) = %g",BiSect(1.55,2,0.1,0.1))
print format$("   BiSect(1.5,  2, .001, .1  ) = %g",BiSect(1.5,2,0.001,0.1))
Waitinput
END
'}-------------------------------------------------------------------------
'Function F(x) HIER EINPORGRAMMIEREN!

proc F :parameters x!

    return cos(x!)

endproc

'-------------------------------------------------------------------------
'Bisec Routine:

proc BISEC :parameters y1!,y2!,epsilon!,epsilon1!

    declare x!,f!,f1!,i&,j&,k&,FLSXT!,ii&
    FLSXT!=-999999999'bedeutet ERROR, keine Nullstelle in diesem Bereich
    Bisec:
    i&=1:j&=1:k&=1:x!=y2!:ii&=0
    Alpha:
    inc ii&':: print "(";ii&;")";
    f!=F(x!) : case abs(f!)<=epsilon!:return x!

    if i&=1:goto "First_val" : elseif i&=2:goto "Succ_val" :endif

        First_val:
        i&=2 : f1!=f! : x!=y1!
        goto "Alpha"
        Succ_val:

        if (f!*f1!)>=0

            if j&=1: Print "\nKeine Nullstelle in diesem Bereich! Error\n";:return FLSXT!

            elseif j&=2: goto "Reg_delta"

            endif

        endif

        if k&=1:goto "Sec_val" : elseif k&=2:goto "Reg_etha" :endif

            Sec_val:
            j&=2:k&=2
            Midpoint:
            x!=(y1!+y2!)/2
            goto "Alpha"
            Reg_delta:
            y2!=x!
            Precision:
            case (abs(y1!-y2!)>=epsilon1!):goto "Midpoint"
            return x!
            Reg_etha:
            y1!=x!
            goto "Precision"

        EndProc

        '-------------------------------------------------------------------------
        'Bisec Routine mit NMAX:

        Proc BiSect

            Declare float y1,y2,epsilon,epsilon1, long Durchlaeufe

            If %PCount = 5

                Parameters float y1_5,y2_5,e_5,e1_5, long NMAX
                y1 = y1_5 : y2 = y2_5 : epsilon = e_5 : epsilon1 = e1_5 : Durchlaeufe = NMAX

            Else

                Parameters float y1_4,y2_4,e_4,e1_4
                y1 = y1_4 : y2 = y2_4 : epsilon = e_4 : epsilon1 = e1_4 : Durchlaeufe = 100

            EndIf

            Declare float x,f0,f1
            Var long Iterationen = 0
            Var float KEINENULL = -999999999.0' bedeutet FEHLER; keine Nullstelle in diesem Bereich
            Var int First = 1
            x = y2

            Repeat

                Inc Iterationen':: print "(";Iterationen;")";
                f0 = F(x)' hier wird die Funktion aufgerufen
                Case Abs(f0) <= epsilon : Return x

                If Iterationen > 1

                    If (f0 * f1) >= 0

                        If First

                            Print "\nKeine Nullstelle in diesem Bereich! FEHLER"
                            Return KEINENULL

                        EndIf

                        y2 = x
                        Case (Abs(y1 - y2) < epsilon1): Return x

                    Else

                        IfNot First

                            y1 = x
                            Case (Abs(y1 - y2) < epsilon1): Return x

                        Else

                            Dec First

                        EndIf

                    EndIf

                    Case (Iterationen < Durchlaeufe) : x = (y1 + y2) / 2

                Else

                    ' der erste Durchlauf
                    f1 = f0
                    x = y1

                EndIf

            Until Iterationen >= Durchlaeufe

            Print "\nZu viele Durchläufe! (";Durchlaeufe;" sind erlaubt)"
            Return x

        EndProc

Computer: Gerät, daß es in Mikrosekunden erlaubt, 50.000 Fehler zu machen, zB 'daß' statt 'das'...

16.05.2021 ▲

Vollständiges Zitat

Absatzweise zitieren Zeilenweise zitieren

p.specht

Andere Umsetzung, von noch früher:

Kompilieren Markieren Separieren

Windowtitle "Bisection-Algorithmus"
'Q: Fortran-IV Source aus 1962, nach XProfan portiert 2015-05 by PS
Cls:Font 2:randomize:AppendMenubar 200,\
"Eingabeln einer Nullstelle per Unterschiedliche_Vorzeichen-Regel"

Proc Fn :Parameters fn&,x!:Declare y!

    Select fn&

        caseof 0: y!=sqr(x!)+exp(x!)-2

        caseof 1: y!=cos(x!)

        caseof 2: y!=3-sin(x!)

        caseof 3: y!=40-exp(x!)

        caseof 4: y!=10

        Otherwise : y!=rnd(1000)
        Endselect :Return y!

    endproc

    Proc ShowFn :Parameters fn& :declare y$

        Select fn&

            caseof 0: y$="y=sqr(x)+exp(x)-2"

            caseof 1: y$="y=cos(x)"

            caseof 2: y$="y=3-sin(x)"

            caseof 3: y$="y=40-exp(x)"

            caseof 4: y$="y=10"

            Otherwise : y$="y=Undefined "
            Endselect :return y$

        endproc

        Proc Bisect

            ' Den Vorzeichenwechsel eingabeln
            Parameters fn&,xl!,xr!,epsx!,epsy!
            Declare yl!,ym!,yr!,xm!
            yl!=Fn(fn&,xl!)

            Repeat

                xm!=(xl!+xr!)*0.5
                ym!=Fn(fn&,xm!)
                case (abs(ym!)-epsy!)<=0:BREAK

                if (yl!*ym!)<0: xr!=xm!

                elseif (yl!*ym!)=0:BREAK

                    else : xl!=xm!

                endif

            Until 0

            Return xm!

        endproc

        '---------------------------------------------------------------
        'Hauptprogramm
        declare fn&,x1!,y1!,x2!,y2!,epsx!,epsy!,fakt!
        '---------------------------------------------------------------
        ' Erfolgs- und performanzkritische Parameter!:
        fn&=0'... Testfunktionsnummer
        epsx!=val("1e-16")
        epsy!=val("1e-16")'...Suchgenauigkeit
        x1!=0 : x2!=0.1'... Startintervall
        fakt!= -1.2'... klappt um und vergrößert Suchbereich
        '---------------------------------------------------------------
        'Eingabeln:
        y1!=Fn(fn&,x1!)
        H1:
        y2!=Fn(fn&,x2!)
        case (y1!*y2!)<0:goto "H4"
        case (y1!*y2!)=0:goto "H5"
        '>0:' gleiche Vorzeichen ==> Bereich erweitern durch ...
        H2:
        x2!=x2!*fakt!'Vergrößerung und Spiegelung an 0-Punkt
        case abs(x2!)<val("1.1e35"):goto "H1"
        H3:
        Print " Sorry, keine Nullstelle gefunden!":goto "H6"
        H4:
        x2!=Bisect(fn&,x1!,x2!,epsx!,epsy!)
        Print "\n Nullstelle der Funktion  ";upper$(ShowFn(Fn&));"\n"
        Print " gefunden bei x = ";format$("%g",x2!);"\n"
        Print " Nullprobe:   y = ";format$("%g",Fn(fn&,x2!));"\n"
        Print " ----------------------------------------------------------"
        H6:
        waitinput
        END