Funktion wird plötzlich nicht geplottet, Integral wird plötzlich kleiner (alle Flächen positiv)

[Mr.Green]

Hey ich habe ein Problem... naja sagen wir mehrerer und ich komme nicht drauf was ich falsch mache. Ich versuche mir grad Python auf dem harten Weg anzueignen, indem ich eine Übung auf der Uni besuche in welcher Beispiele programmiert werden müssen, jedoch keine Hilfe zu den Sprachen gegeben werden kann, da der Prof Python auch nicht wirklich beherrscht, dementsprechend verzweifelt bin ich langsam, da ich seit mehreren Tagen nicht zu einem Ergebnis komme ^^

Nun zu den Problemen: Plotte ich meine Funktion auf einem Intervall von 1 bis 3 klappt alles (stelle ich über Linspace ein) Plotte ich sie jedoch von 0 bis 10 (wie jetzt gerade), verschwindet meine Funktion und ich habe in der Ausgabe nur mehr eine gerade Linie...

Mein Zweites Problem (ich hoffe es hängt mit dem Ersten zusammen) Integriere ich mit integrate.quad meine Funktion mit den Grenzen von 1 bis 3 so bekomme ich einen akzeptablen Wert. Dies ändert sich jedoch wenn ich meine rechte Grenze in Richtung 1000 verschiebe (wie hier jetzt) dann kommt als Ergebnis 0 bzw. ein sehr kleiner Wert heraus. Das kann aber nicht sein, da all meine Werte der Funktion positiv sind (außer ich hab mich ca. 10 Mal beim durchdenken vertan, was ich mittlerweile auch nicht mehr ausschließen möchte)

Ich hoffe ihr könnt mir helfen und mich vor meinem persönlichen "schwarzen" Freitag bewahren.
LG Christoph

Code: Alles auswählen

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
x = np.linspace(0,10,10000)
def p(l,n,S):
        return (n*l)**S*np.exp(-n*l)*1/(1+0.25*l)
def lnp(l,n,S,K):
    return np.exp(S*np.log(n*l)-n*l-np.log(1+0.25*l)-K)
n=60
S=124
K = np.amax(np.log(p(x,n,S)))
print(K)
def f(a,n,S):
        return p(a,n,S)*np.exp(-K)
y =lnp(x,n,S,K)
plt.plot(x,y)
plt.show()
b=integrate.quad(lnp,1,1000,args=(n,S,K))
print(b)

Kurz zur Erklärung. Die Funktion p ist meine Ausgangsfunktion, welche jedoch teils sehr große Werte annehmen kann, daher habe ich den ln genommen und mein Maximum gesucht welches ich als K ausgebe. Mit diesem K senke ich mein Maximum auf 0 ab, damit ich, wenn ich wieder die Exponentialfunktion dieser ln-Funktion nehme meinen Maximal-Wert bei 1 habe. Danach plotte ich diese Funktion und bilde das Integral darüber. (möchte eigentlich von 0 bis unendlich integrieren aber.... es geht halt nicht einmal bis 1000 :/
LG Christoph

[BlackJack]

@Mr.Green: Schau Dir doch einfach mal die Zwischenergebnisse an. `K` beispielsweise. Das hat den Wert ∞, was wahrscheinlich nicht das ist was Du willst‽

[Mr.Green]

K sollte jedoch nie unendlich werden. Plotte ich nur die Funktion aus welcher ich dann das Maximum suche so erreicht K nur einen Maximalwert.

Code: Alles auswählen

K = np.amax(np.log(p(x,n,S)))
print(K)
R = np.log(p(x,n,S))
plt.plot(x,R)
plt.show()

K ist auch normal solange ich mich bei x=np.linspace unterhalb von 5 befinde, erst danach wird K plötzlich unendlich.

[Sirius3]

@Mr.Green: Du rechnest mit Zahlen, die nur eine endliche Genauigkeit haben. Irgendwas hoch 124 sieht dann ziemlich schnell wie ∞ aus.

[Mr.Green]

stimmt!! genau deswegen nehme ich ja später den ln der Funktion, habe aber vollkommen übersehen, dass mein K ja noch mit der ursprünglichen Funktion berechnet wird...Jetzt geht's ^^... ich war ja echt blind
Danke
LG Christoph

[Mr.Green]

muss korrigieren... es geht beinahe... meine Integralgrenze kann ich immer noch nicht höher als 1000 setzen...

[BlackJack]

@Mr.Green: Du nimmst nicht den Logarithmus der Funktion sondern den Logarithmus vom Ergebnis der Funktion. Wenn die Funktion also irgendwelche Unendlich-Werte liefert, dann bekommst Du die mit dem Logarithmus nicht wieder zurück, denn dann ist es zu spät.

Die Grenzen der Gleitkommatypen kann man sich übrigens mit `numpy.finfo()` geben lassen:

Code: Alles auswählen

In [16]: np.finfo(float)
Out[16]: finfo(resolution=1e-15, min=-1.7976931348623157e+308, max=1.7976931348623157e+308, dtype=float64)

In [17]: np.finfo(np.float96)
Out[17]: finfo(resolution=1e-18, min=-1.18973149536e+4932, max=1.18973149536e+4932, dtype=float96)

Edit: Ob es den letzten Typ gibt, hängt von der Plattform ab auf der es läuft. 64-Bit-Floats gibt es auf jeder Plattform weil das C's ``double``-Datentyp entspricht und die notfalls in Software gerechnet werden wenn die Hardware die nicht ”nativ” bietet.

[Mr.Green]

Für mein K habe ich bis jetzt erst vom Endergebnis den log genommen, dadurch hat es mir für K sehr schnell unendlich ausgegeben, jedoch für mein Integral benutze ich immer schon die log-Funktion, in diesem Fall dürften die Werte also nicht zu groß werden.
(Wäre es vielleicht allgemein besser wenn ich das Integral über sympy löse?)
Danke für den Tipp, bin ja wirklich sehr neu bei Python und freue mich immer wenn ich dann parallel zu Tutorien noch so etwas aufschnappen kann was mir hilft.

[Sirius3]

@Mr.Green: das kommt halt darauf an, was Du machen willst, und was Du lernen willst. Wenn es um die numerische Auswertung von Funktionen geht, dann sollst Du ja wahrscheinlich über die Stolpersteine fallen und Dir Gedanken machen, wann eine numerische Integration realistische Werte liefert, wann nicht und was man dagegen tun kann. Wie man erkennt, dass etwas falsch läuft, wie man das verhindert, usw. Falls das nicht das Ziel ist, macht es hier mehr Sinn, die Funktion händisch zu integrieren und dann zu prüfen, ob sympy auch auf das Ergebnis kommt.