Fortgesetzte Bisektion bei bestimmter Funktion

Schaut euch mal die folgende Funktion an:

Figure_1

Wie nicht schwer zu erkennen ist, ist irgendwo zwischen 0 und 2 ein Minimum. Wenn man die Funktion ableiten würde, dann wäre das Minimum schnell erkennbar.

Aber kann man das Minimum der Funktion auch mit dem Intervallhalbierungsverfahren bestimmen? Wenn ja, wie?

Die Funktion (f4) ist folgendermaßen definiert:

import math

# importing the required module 
import matplotlib.pyplot as plt

def f1(h:float, v:float):
    return math.sqrt(v / (math.pi * h))

def f2(h:float, r:float, k1:float, k2:float):
    return 2 * math.pi * r * h * k1 + 2 * math.pi * r * r * k2;

def f3(h:float, v:float, k1:float, k2:float):
    return f2(h, f1(h, v), k1, k2)

def f4(h:float):
    v:float = 50
    k1:float = 0.05
    k2:float = 0.01
    return f3(h, v, k1, k2)

def plot_f4():
    x:[float] = []
    y:[float] = []
    for i in range(1, 160):
        h:float = i / 10.0
        x.append(h)
        y.append(f4(h))
    # plotting the points
    plt.plot(x, y)
    # naming the x axis
    plt.xlabel('x - axis')
    # naming the y axis
    plt.ylabel('y - axis')
    # giving a title to my graph
    plt.title('My first graph!')
    # function to show the plot
    plt.show()

plot_f4()

Danke für jeden Hinweis!

Edit: Das ist - btw. - nicht ganz aus der Luft gegriffen… Mit Hilfe der Funktion kann das optimale Verhältnis zwischen Höhe und Durchmesser einer Getränkedose in der Verpackungsherstellung berechnet werden.

Für Funktionen mit einem einzigen Minimum im Ausgangsbereich sollte folgendes funktionieren:

Gegeben sei das Intervall a,b. Dann berechnet man nicht den Mittelpunkt, sondern drittelt das Intervall, also c = 2a/3 + b/3 und d = a/3 + 2b/3.

  • wenn f( c ) < f( d ), mache mit dem Intervall a,d weiter
  • wenn f( d ) < f( c ), mache mit dem Intervall c,b weiter
  • wenn f( c ) = f( d ), mache mit dem Intervall c,d weiter

Danke Landei!

Hier noch zwei Beispiele:
Es sollen 50ml Flüssigkeit verpackt werden. Ein Quadratzentimeter der Mantelfläche/Seiten kostet 1 Cent. Ein Quadratzentimeter des Bodens und des Deckels kosten 0,5 Cent. Wie ist die Höhe und der Radius zu wählen, damit die Gesamtkosten am kleinsten sind? (Hier ist das Verhältnis 1:1 gut.)

Kosten qcm Boden und Deckel der Dose in Euro:
0.005
Kosten qcm Mantelfläche der Dose in Euro:
0.01
Benötigtes Volumen der Dose in Kubikzentimeter (Milliliter):
50
Optimale Höhe         : 2.5250000000000923
Optimale Breite       : 5.021221390476676
Optimaler Radius      : 2.510610695238338
Optimales Volumen     : 50.00000000000001
Optimale Gesamtfläche : 79.43490678277249
Optimale Gesamtkosten : 0.596329265847534
Gesamtkosten (h-0.5)  : 0.6036130370314515
Gesamtkosten (h+0.5)  : 0.6012552589350098

Es sollen 50ml Flüssigkeit verpackt werden. Ein Quadratzentimeter der Mantelfläche/Seiten kostet 1 Cent. Ein Quadratzentimeter des Bodens und des Deckels kosten 2 Cent. Wie ist die Höhe und der Radius zu wählen, damit die Gesamtkosten am kleinsten sind? (Hier ist das Verhältnis 1:4,01 gut.)

Kosten qcm Boden und Deckel der Dose in Euro:
0.02
Kosten qcm Mantelfläche der Dose in Euro:
0.01
Benötigtes Volumen der Dose in Kubikzentimeter (Milliliter):
50
Optimale Höhe         : 6.349999999999993
Optimale Breite       : 3.1663082961488658
Optimaler Radius      : 1.5831541480744329
Optimales Volumen     : 49.99999999999999
Optimale Gesamtfläche : 78.91307459791669
Optimale Gesamtkosten : 0.9466110609397969
Gesamtkosten (h-0.5)  : 0.9481528387577821
Gesamtkosten (h+0.5)  : 0.9480182048212791

Landei, was hältst du von dieser Methode? Ist das Intervall-Halbierung?

def calc_optimum_hoehe():
    v:float = 50
    k1:float = 0.01
    k2:float = 0.005
    e:float = 0.1
    h:float = 0.1
    min:float = 100000
    for i in range(500):
        fy:float = f3(h, v, k1, k2)
        if fy <= min:
            min = fy 
            h += e 
        else:
            e /= 2.0
            h -= e 
    return h 

print(calc_optimum_hoehe())

Könnte funktionieren, ist aber meiner Meinung nach keine Intervall-Halbierung. Intervall-Halbierung ist: Intervall halbieren, schauen, und dann mit der „richtigen“ Hälfte weitermachen.

Wenn du unbedingt Invervall-Halbierung haben willst: Intervall halbieren, Funktionswerte an den Mittelpunkten der Teilintervalle vergleichen, mit dem Teilintervall weitermachen, das den kleineren Wert hat.

Das ist jetzt doch ein Form von Simulated Annealing oder liege ich jetzt ganz falsch? :clown_face:

Beispiel:
https://de.wikipedia.org/wiki/Simulated_Annealing#/media/Datei:SimAnnealingLandschaft.png

Nein, Simulated Annealing ist komplexer, da steigt die lokale Temperatur auch manchmal, um nicht in lokalen Minima steckenzubleiben.

Simulated Annealing ist kein geschützter Begriff. Du kannst das nennen, wie du willst. Man kann auch „Kurvenfitten“ (oder sogar „einen Haufen if-Abfragen“) als „Künstliche Intelligenz“ bezeichnen, und viel Geld damit verdienen, und niemand kann einem was.

Ansonsten eine Antwort, die dir vom Stil her zusagen könnte:

public class BisectionTest
{
    public static void main(String[] args)
    {
        DoubleUnaryOperator î = x -> x * x * x + 3 * x * x - 6 * x - 8;

        double í = i(-3.5, 1.2, 1e-4, î);
        System.out.println(í);
    }

    private static double i(double í, double ì, double i, DoubleUnaryOperator î)
    {
        return Math.abs(í-ì)<i?(í+ì)/2:î.applyAsDouble(í)<î.applyAsDouble((í+
            ì)/2)?i(í,(í+ì)/2,i,î):i((í+ì)/2,ì,i,î);
    }
}

Ich weiß nicht, ob sie immer die richtige Lösung liefert, aber lies’s sie dir mal durch, und schau, ob du einen Fehler findest.

Was soll das? Ich fühle mich von dir veräppelt. :wink:

Spricht etwas dagegen, nicht a, b, c, d und e zu verwenden?

Sorry, ich weiß, was Marco13 meint:

import math
import sys

# importing the required module 
import matplotlib.pyplot as plt

def gesamtflaeche(hoehe:float, radius:float):
    return 2 * math.pi * radius * (radius + hoehe)

def volumen(hoehe:float, radius:float):
    return math.pi * radius * radius * hoehe

def radius(hoehe:float, volumen:float):
    return math.sqrt(volumen / (math.pi * hoehe))

def gesamtkosten_1(h:float, r:float, kosten_mantel:float, kosten_boden_deckel:float):
    return 2 * math.pi * r * h * kosten_mantel + 2 * math.pi * r * r * kosten_boden_deckel;

def gesamtkosten_2(h:float, v:float, kosten_mantel:float, kosten_boden_deckel:float):
    return gesamtkosten_1(h, radius(h, v), kosten_mantel, kosten_boden_deckel)

def gesamtkosten_3(hoehe:float):
    volumen:float = 500
    kosten_mantel:float = 0.01
    kosten_boden_deckel:float = 0.005
    return gesamtkosten_2(hoehe, volumen, kosten_mantel, kosten_boden_deckel)

def plot_gesamtkosten_3():
    x:[float] = []
    y:[float] = []
    for i in range(30, 70):
        hoehe:float = i / 10.0
        x.append(hoehe)
        y.append(gesamtkosten_3(hoehe))
    # plotting the points
    plt.plot(x, y)
    # naming the x axis
    plt.xlabel('Hoehe')
    # naming the y axis
    plt.ylabel('Kosten')
    # giving a title to my graph
    plt.title('Kosten vs. Hoehe!')
    # function to show the plot
    plt.show()

def calc_optimum_hoehe():
    volumen:float = 500
    kosten_mantel:float = 0.01
    kosten_boden_deckel:float = 0.005
    e:float = 0.1
    hoehe:float = 0.1
    min:float = sys.float_info.max
    for i in range(500):
        fy:float = gesamtkosten_2(hoehe, volumen, kosten_mantel, kosten_boden_deckel)
        if fy <= min:
            min = fy 
            hoehe += e 
        else:
            e /= 2.0
            hoehe -= e 
    return hoehe 

print(calc_optimum_hoehe())
print(radius(calc_optimum_hoehe(), volumen=500))
print(volumen(calc_optimum_hoehe(), radius(calc_optimum_hoehe(), volumen=500)))
print(gesamtflaeche(calc_optimum_hoehe(), radius(calc_optimum_hoehe(), volumen=500)))

plot_gesamtkosten_3()

… besser?

Die Bisektion ist eher für die Bestimmung einer Nullstelle gedacht. Zur Bestimmung des Minimums ist sie eher ungeeignet. Das klappt nicht immer.

Danke für Deine Antwort! … Das ist mir auch schon aufgefallen, daher mein modifiziertes „Simulated Annealing“, wobei sich die Temperatur stets abkühlt.

Bestimme einfach die Nullstelle der Ableitung. (Hier könnte ein :clown_face: stehen, aber das steht hier ganz bewußt nicht…)

Ja, es ist die Nullstelle der Ableitung, aber diese wollte ich ja bewusst nicht herleiten… Sondern, ich habe nach einem „allgemeinen“ Verfahren gefragt, wenn ein Minimum oder Maximum gesucht werden soll. Hast du 'ne Idee?

Edit: Die Nullstelle der Ableitung könnte sogar rechnerisch „exakt“ berechnet werden, imo.

Der Einwand von nolle zielte vemutlich (!?) darauf ab, dass man bei der Suche auch in ein lokales Extremum rutschen kann - falls das gemeint war, ergibt es aber nur bedingt Sinn, weil das gleichbedeutend damit ist, dass die Ableitung mehrere Nullstellen hat, und man die „falsche“ erwischt.

Die Ableitung zu bilden (algebraisch oder numerisch) und dort die Nullstelle zu suchen ist ein allgemeines Verfahren. Da du sowieso irgendwie nach Steigungen schauen musst, arbeitest du im Prinzip „implizit“ mit der Ableitung, egal wie sehr du es vermeiden willst.

„Simulated Annealing“ schaut nicht nach der Steigung/Ableitung, das ist ja das Schöne. :clown_face:

Zitat:
„Simulated annealing (SA) is a probabilistic technique for approximating the global optimum of a given function.“

probabilistic = wir probieren etwas aus. :slight_smile: Man sucht sich im Prinzip systematisch zufällige Nachbarn aus und schaut, ob man tiefer als bisher ist - meist hat man Glück und landet beim globalen Minimum…

Da hast du deine Steigung.

Tada! Simulated Annealing:

def sa_gesamtkosten_3():
    alpha = 0.1
    initial_temp = 100
    final_temp = 1
    current_temp = initial_temp
    solution = current_temp / 2.0
    while current_temp > final_temp:
        new_h = random.gauss(solution, current_temp / 2.0)
        while new_h <= 0:
            new_h = random.gauss(solution, current_temp / 2.0)
        costs_x = gesamtkosten_3(solution)
        costs_y = gesamtkosten_3(new_h)
        if costs_y < costs_x:
            solution = new_h
        elif random.random() > math.exp((costs_y - costs_x) / -initial_temp):
            solution = new_h
        current_temp -= alpha
        print(current_temp, new_h, solution)
    return solution

Ich bin mir allerdings bei der Wahl des Nachbarn (random.gauss(solution, current_temp / 2.0)) und beim Mogeln (if random.random() > math.exp((costs_y - costs_x) / -initial_temp)) (sollte auch mal eine schlechtere Lösung übernommen werden?) unsicher, ob die Terme/Variablen richtig sind.

Wenigstens spuckt es ein richtiges Ergebnis aus. :slight_smile: