Numerische Integration

Wenn du dir nicht sicher bist, in welchem der anderen Foren du die Frage stellen sollst, dann bist du hier im Forum für allgemeine Fragen sicher richtig.
Heini
User
Beiträge: 21
Registriert: Donnerstag 12. April 2007, 10:48

Hallo.
Ich habe derzeit ein sehr unfreundliches 3dimensionales Integral numerisch auszurechnen. Die tplquad-Funktion von scipy ist dafür leider zu langsam oder besser gesagt, sie liefert kein Ergebnis. Meinem Chef gelingt es jedoch das Integral in Fortran auszurechnen, indem er eine Integration mit Hilfe der Gaußintegration bestimmter Ordnung ausführt und die anderen beiden indem er die Simpsonmethode geschickt darauf los lässt. Das gleiche habe ich dann in Python probiert, wobei ich für die Simpsonmethode den Code auf (die gesamte Zeile ist die Adresse)

http://en.wikipedia.org/wiki/Adaptive_Simpson's_method

benutze. Leider weiß ich nicht wie man die dort beschriebene Routine effektiver gestallten kann und ob es überhaupt daran liegt. Am Ende gelingt es mir dann auch nicht irgendein Ergebnis zu erhalten und das obwohl ich noch nicht einmal die vollständige Formel benutze.
Die Profis unter euch werden schon gemerkt haben: "...der hat keine Ahnung vom Programmieren." Sehr richtig. Dennoch muss das Mistviech irgendwie berechnet werden. Wenn also jemand eine unglaublich schnelle und genaue Methode parat hat wie man solche unliebsamen Weggefährten summieren kann, dann wäre ich für einen guten Rat sehr dankbar.
CM
User
Beiträge: 2464
Registriert: Sonntag 29. August 2004, 19:47
Kontaktdaten:

Hoi & willkommen im Forum,

davon ab, daß die Simpsonsche-Regel ebenfalls in scipy enthalten ist: Hier wird Dir niemand helfen können, wenn Du uns nicht die Form Deiner Daten zeigst. (Handelt es sich um eine 3D-Matrix in numpy?) Außerdem wäre es nicht schlecht Deinen bisherigen Versuch zu sehen oder vielleicht sogar den Ansatz Deines Chefs. Und wenn Du möchtest, daß mir mehr als Raten brauchen wir tatsächlich die Daten zum Testen - natürlich nur, wenn es nicht zu groß ist (und am Besten noch das Ergebnis ;-) ).

Gruß,
Christian
Heini
User
Beiträge: 21
Registriert: Donnerstag 12. April 2007, 10:48

Ich habe bereits die Gegenbauerpolynome per Hand eingefügt um den Rechenaufwand zu verringern. Die Integrationsgrenzen sind:
cos(chi)=-1..+1, k=0..100, k4=-5..+5; die Grenzen für k und k4 müssen nicht unbedingt größer gewählt werden, aber auch nicht kleiner.
Wenn man das einfach so eingibt, und tplquad auf die vollen Intervalle loslässt, hat man glaube ich keine Chance.
Die cos(chi)-integration kann man sehr gut mit der Gauß-Methode mit n=96 ausführen, in dieser Variable ist die Funktion "gutartig". Aber auch wenn man danach dblquad darauf ansetzt bekommt man nicht unbedingt schnell eine Lösung.
Zuletzt geändert von Heini am Donnerstag 12. April 2007, 14:02, insgesamt 1-mal geändert.
CM
User
Beiträge: 2464
Registriert: Sonntag 29. August 2004, 19:47
Kontaktdaten:

Heini hat geschrieben:Also, ich poste das mal wie folgt. Zuerst die Formel(n) mit den Konstanten, danach meine beiden Versuche ohne die Formel. Zugegeben, beide Programme sehen aus wie Kraut und Rüben, aber naja...
In der Tat:
http://www.python-forum.de/faq.php#21

habe gerade einen Termin - schau nachher nochmal drauf.

Christian
Heini
User
Beiträge: 21
Registriert: Donnerstag 12. April 2007, 10:48

kleiner Fehler:
def g(k4,k): -> def g(k4,k,i):
(es ist nur ca. 1/32 der Gesamtformel, daher stehen da noch "überflüssige Indizes bei g() und Z())
Heini
User
Beiträge: 21
Registriert: Donnerstag 12. April 2007, 10:48

ah, sorry...na was solls, ich pack einfach mal einen Versuch ganz rein. Vielleicht ist ja ein Admin so nett den obigen Salat zu entfernen.

Code: Alles auswählen


""" This program calculates the electromagnetic formfactor for the scalar
    1S0-channel in static approximation.
"""

from scipy.integrate import tplquad, dblquad
from numpy import *
from scipy import sqrt
from scipy.special.orthogonal import p_roots

# the function -------------------------------------
def g(k4, k, i):
    x = sqrt(k4**2+k**2)                            # =abs(k_E)
    t = k4/x                                        # =cos(theta)
    if i == 1:
        g1 = 0.17926/(0.42449**2+x**2)
        g2 = -0.01057*x**2/(0.28844**2+x**2)**2*(4*t**2-1)
        g3 = 0.00122*x**4/(0.37345**2+x**2)**4*(16*t**4-12*t**2+1)
        g4 = -0.00025*x**5/(0.41686**2+x**2)**5*(64*t**6-80*t**4+24*t**2-1)
    elif i==2:
        g1 = 0.14667/(0.39778**2+x**2)
        g2 = -0.0085*x**2/(0.27192**2+x**2)**2*(4*t**2-1)
        g3 = 0.0009*x**4/(0.35564**2+x**2)**4*(16*t**4-12*t**2+1)
        g4 = -0.00015*x**6/(0.4008**2+x**2)**6*(64*t**6-80*t**4+24*t**2-1)
    elif i == 3:
        g1 = 0.04037*x/(0.5978**2+x**2)**2
        g2 = -0.00298*x**3/(0.42496**2+x**2)**3*(12*t**2-2)
        g3 = 0.00056*x**4/(0.47242**2+x**2)**4*(80*t**4-48*t**2+3)
        g4 = -0.00013*x**5/(0.49936**2+x**2)**5*(448*t**6-480*t**4+120*t**2-4)
    elif i == 4:
        g1 = -0.0033*x**2/(0.00051**2+x**2)**2*(4*t)
        g2 = 0.00051*x**3/(0.66046**2+x**2)**3*(32*t**3-12*t)
        g3 = -0.00011*x**4/(0.65207**2+x**2)**4*(192*t**5-160*t**3+24*t)
        g4 = 0.00003*x**5/(0.64625**2+x**2)**5*(1024*t**7-1344*t**5+480*t**3-40*t)
    return sqrt(2/pi)*(g1+g2+g3+g4)

def Z(k4, k, cos_theta, i, j):
    t = cos_theta
    k0 = 1j*k4
    if i == 1:
        if j == 1:
            res = .5*M**3*k0+2*k*M**2+k**2*M**2-M**3*alpha*(k*t+M*alpha*.5)*.25\
                -2*k0**3*M*etha-k0**2*M**2*etha+3*m**2*M**2*etha+2*k**2*M*k0\
                +2*m**2*M*k0+M**3*alpha*k*.5-k0**2*M**2-2*k*t**2*M**2\
                +M**4*etha*.25+2*k*t*(k*t+M*alpha*.5)*M**2*etha\
                +k0**2*M*alpha*(k*t+M*alpha*.5)-m**2*M*alpha*(k*t+M*alpha*.5)\
                -k**2*M*alpha*(k*t+M*alpha*.5)-M**2*k0*alpha*k\
                +2*k**2*M*k0*etha-2*k*t**2*M**2*etha+2*m**2*M*k0*etha\
                +2*k*t*(k*t+M*alpha*.5)*M**2-2*k0**3*M+2*k*M**2*etha\
                +M**3*k0*etha*.5+M**4*.25+3*m**2*M**2+k**2*M**2*etha
    return res

def F(k, k4, cos_theta):
    """ The function to be integrated.
    """
    epsk = sqrt(k**2+m**2)
    k01k02 = M**2*.25*(1+4*etha)-(M**2*alpha**2+2*M*alpha*k*cos_theta+epsk**2)
    prop1 = 1.0/((M*.5+epsk)**2+k4**2)
    prop2 = 1.0/((M*.5-epsk)**2+k4**2)
    prop3 = (-k4**2+k01k02-1j*k4*M*(1+4*etha))/((-k4**2+k01k02)**2\
        +k4**2*M**2*(1+4*etha)**2)
    res = g(k4, sqrt(k**2+q**2*.25+k*q*cos_theta), 1)*g(k4, k, 1)
    res *= k**2*Z(k4, k, cos_theta, 1, 1)*prop1*prop2*prop3/(2*pi)**3
    return res
#-------------------------------------------------
def G(k, k4):
    """ Performing the inner integration over cos(theta) by gaussian method
        for fixed k, k4.
        This function is used in combination with dblquad instead of tplquad.
        F(k, k4, cos_theta) -> G(k, k4)
    """
    return sum(w*F(k, k4, x))
#-------------------------------------------------

# input parameters
m = 1.0                                             # mass
M = 1.974490569493619                               # mass

x,w = real(p_roots(96))
for i in xrange(21):                              # vary Q
    # more constants
    Q = sqrt(0.2)
    etha = Q**2/(4*M**2)
    alpha = 2*sqrt(etha*(1+etha))
    q = M*alpha                                   # abs(vec{q})
    # integration is performed on 2 intervals
    # 1) integration using dblquad and gaussian method
    ff = dblquad(G, -5.0, 5.0, lambda x: 0.0, lambda x: 3.0)[0]
    ff += dblquad(G, -5.0, 5.0, lambda x: 3.0, lambda x: 100.0)[0]
##    # 2) integration using tplquad
##    ff = tplquad(F, -1.0, 1.0, lambda z: -5.0, lambda z: 5.0, lambda x,y: 0.0, lambda x,y: 3.0)[0]
##    ff += tplquad(F, -1.0, 1.0, lambda z: -5.0, lambda z: 5.0, lambda x,y: 3.0, lambda x,y: 100.0)[0]
    print 'Q^2=%.3f  ff=%.18f '%(Q**2, ff)

F,Z und g sind im Grunde die besagte Formel. Der K(r)ampf mit der Integration beginnt danach. Hier habe ich bereits versucht sowohl die Intervalle vorher zu unterteilen, damit dblquad nicht erst von "ganz oben" anfangen muss, als auch die cos(chi)-Integration mit "manueller"-Gaussintegration durchzuführen.
Das Integral muss nicht unbedingt in wenigen Sekunden gelöst werden, zumal es ja mehrmals ausgerechnet werden soll, aber endlich sollte die Zeit schon bleiben. Wenn für ein bestimmtes Q das Integral nach ein paar Minuten berechnet werden würde, wäre ich schon sehr zufrieden.
Zuletzt geändert von Heini am Freitag 13. April 2007, 07:09, insgesamt 3-mal geändert.
Benutzeravatar
mq
User
Beiträge: 124
Registriert: Samstag 1. Januar 2005, 19:14

Erwartest du ernsthaft, dass jemand diesen Code versteht und darin auch noch Bugs findet?
Zuletzt geändert von mq am Donnerstag 12. April 2007, 15:43, insgesamt 1-mal geändert.
lunar

lumax hat geschrieben:Erwartest du, dass jemand diesen Code versteht und darin auch noch Bugs findet?
Na, offenbar schon, sonst hätte er ihn ja nicht gepostet...
CM
User
Beiträge: 2464
Registriert: Sonntag 29. August 2004, 19:47
Kontaktdaten:

Hoi,

puh, glaube daran werde ich mir die Zähne ausbeißen. ;-) Vielleicht magst Du unter der scipy-Mailingliste nochmal nachfragen?

Aber ich habe durchaus noch weitere Fragen/Anmerkungen, die helfen könnten das Problem einzugrenzen, bzw. den Code etwas besser zu gestalten (in willkürlicher Reihenfolge):
- Wozu brauchst Du die Doppelschleife? Du greifst weder in dblquad noch in tplquad auf n (oder xk) zu und Du summierst auch nichts auf oder so.
- Ohne entsprechende Docstrings und Kommentare wird kaum jemand Deinen Code verstehen können - Du selber in ein paar Wochen auch nicht mehr. Dann wird evtuelles Debuggen richtig spaßig.
- Es gibt in Python keinen Grund alle Parameter im Kopf zu deklarieren und das auch noch ohne Kommentar wofür diese Parameter stehen. Besser ist es - insbesondere damit Du später Deinen Code verstehst - die Paramter dort zu definieren wo sie gebraucht werden und sie mit einem Kommentar zu versehen was sie eigentlich sind.
- Der Code ließt sich wie eine Mischung aus Fortran und Python - kann das sein? ;-) Jedenfalls nicht sehr übersichtlich.
- Die Toleranzgrenzen finde ich extrem klein. Reichen die Defaultwerte nicht? (im nächsten Punkt nehme ich den Defaultfehler und der Fehler ist kleiner als Deine angegebene Toleranz!)
-

Code: Alles auswählen

 ff = tplquad(F,-5.0, 5.0, mink,maxk, lambda x,y: -1.0, lambda x,y: 1.0)
            print ff
gibt mir für n==0: (3.6465208203935958e-06, 7.8470743910729225e-11) in redlicher Zeit auf meinem alten Laptop. Ist das korrekt?
-

Code: Alles auswählen

print 'Q^2=',Q**2,'ff=%.18f'%ff
würde ich so schreiben

Code: Alles auswählen

print 'Q^2=%.3f  ff=%.18f '% (Q**2, ff)
- das nur nebenbei.
- Wieso läßt Du ff in der Doppelschleife immer wieder fallen und hebst das Ergebnis nicht auf? (Du merkst: Ichhabe das Problem noch nicht nachvollzogen.)
- Deine Formeln sind wahre Ungetüme. In Python geht:

Code: Alles auswählen

x *= a # anstatt: x = x * a
x += a # anstatt: x = x + a
und man sollte ein Leerzeichen nach dem Komma machen - nur für die Übersichtlichkeit, Du wirst sehen es hilft -, aber das sind Marginalien.
- Ich glaube, daß Z() und F() hinsichtlich der Operatorpräferenz richtig sind. Bei g() bin ich mir nicht sicher. Wie ist zum Beispiel

Code: Alles auswählen

g2=-0.01057*x**2/(0.28844**2+x**2)**2*(4*t**2-1)
zu verstehen? Gilt
(x**2)/(0.28844**2+x**2) (was da steht) oder
x**(2/(0.28844**2+x**2))? etc.

So, jetzt bin ich gerade angefangen, muß aber für heute wieder aufhören - die Arbeit ruft.

Gruß,
Christian
BlackJack

In `mink()` und `maxk()` wird wohl das `n` benutzt. Insgesamt ein schönes Beispiel warum man Argumente anstelle von globalen Variablen benutzen sollte.
Heini
User
Beiträge: 21
Registriert: Donnerstag 12. April 2007, 10:48

Vielen Dank erstmal für die vielen hilfreichen Hinweise.
Zu meiner Verteidigung muss ich sagen, dass ich den Code selbst auch fürchterlich unordentlich finde, aber er ist dann nach zahllosem ausprobieren irgendwie zu dem Monster angewachsen das er jetzt ist. Aber nun zu den vielen Punkten.
- Die Doppelschleife habe ich eingeführt, um die Intervalle für k und k4 aufzuteilen bevor ich tplquad oder dblquad darauf loslasse. Ich dachte, dass sie dann vielleicht weniger zu tun hätten. Ich greife dann in den Grenzfunktionen (z.B.: mink4 und maxk4) auf diese integer zurück. Ich wusste einfach nicht wie ich das sonst hinbiegen soll, da diese Grenzfunktionen ja nur eine Variable nehmen sollen. Aber du hast recht, ich lasse alle Intervallintegrale einfach von dannen schwimmen.
- Manche der Parameter brauche ich in mehreren Funktionen und ich wollte die nicht andauernd übergeben müssen (wie beispielsweise die Massen) oder nicht immer neu berechnen (wie x,w=real(p_roots)). Aber ich dachte bei err (Fehler) eigentlich daran dass diese Variable dann nicht jedes mal neu definiert werden muss. Warum macht das in Python keinen Sinn? (eigentlich kenne ich nur Python, daher ist jede Ähnlichkeit mit einem Fortrancode reiner Zufall). Für den späteren Gebrauch ist es auch nicht schlecht wenn ich schnell sehe welche Parameter ich benutzt habe.
- Was hat das mit der lambda-funktion auf sich? Kann ich damit die Intervalle aufteilen, so wie ich es eigentlich vorhatte? Das Ergebnis scheint zu stimmen, für n=0 wäre das aber nur das Intervall [0,0.1].
- (x**2)/(0.28844**2+x**2) das soll da stehen. Die Funktion ist sehr scharf.

Kann man vielleicht etwas Zeit gewinnen, wenn man die Reihenfolge der Integrationen vertauscht? Vielleicht fällt es ihm einfacher, wenn die innerste Integration über cos(theta) ist und die äußerste über k; der größe der Intervalle nach geordnet?

Ich versuche mal alles einzubauen und insbesonder den code übersichtlicher zu gestallten.
Vielen Dank nochmal.
Gruß, Thomas
CM
User
Beiträge: 2464
Registriert: Sonntag 29. August 2004, 19:47
Kontaktdaten:

... so, habe doch noch eine Pause - das Problem läßt mir keine Ruhe ;-).

@lumax und lunar: Der Code ist nicht besonders schön, aber Heini hat ja schon geschrieben Anfänger zu sein. Im Übrigen versucht man in der Numerik oftmals Funktionsaufrufe zu vermeiden, um den Code zu beschleunigen. Das führt u. U. zu wahnwitzig langen Formeln.

@BJ: Ups, danke das habe ich übersehen. Aber solche "Designfehler" tendiere ich zu übersehen, weil sie mir kaum unterlaufen.

Aber das ändert nichts am Problem, daß der Code nicht kommentiert ist und damit kaum zu verstehen ist, warum da eine Doppelschleife durchlaufen wird.

Weitere Anmerkungen:
- Um den Code zu beschleunigen kannst Du auch darüber nachdenken so wenige Divisionen wie möglich durchzuführen. Ihr werdet lachen: Multiplizieren in C und auch Python (und auch Fortran) ist schneller als Dividieren und bei einigen tausend Durchgängen fällt so etwas ins Gewicht - insbesondere bei O(n**2)!
- "return (x)" kann man auch "return x" schreiben. Damit sparst Du wenig Zeit, aber immerhin. Außerdem: Pythons automatisches tuple unpacking ist nicht nur praktisch, es macht den Code auch übersichtlich.
- Wichtiger noch: Greife unbedingt BlackJacks Kommentar auf!

Gruß,
Christian
CM
User
Beiträge: 2464
Registriert: Sonntag 29. August 2004, 19:47
Kontaktdaten:

Ups, zu langsam ...
Warum macht das in Python keinen Sinn?
Was ich meinte ist folgendes:

Code: Alles auswählen

def foo(*args):
   pass

if __name__ == "__main__":
    # Erklärung zum Parameter hier
    p1 = 3.1459...
    # Erklärung zum Parameter hier
    p2 = 2.718...
    
    # Hier erklären was Du tun möchtest
    value = foo(p1, p2)
Ließt sich doch übersichtlicher, oder? Vorteil: Die Paramter stehen direkt über der Funktion, die sie aufruft. Der nächste Schritt für Dich wäre Skripte etwas dynamischer zu machen und Parameter aus einem File oder direkt vom Nutzer zu nehem. Auch hierbei empfiehlt es sich diese Dinge übersichtlich abzuwickeln. Pfeife drauf, daß Du Parameter mehrmals übergeben mußt. Vorteil wenn Du das tust: Du sicherst Dich gegen Deine eigenen Fehler ab. Das erschließt sich Dir, wenn Du länger codest.
Was hat das mit der lambda-funktion auf sich?
Nichts, aber was machen Deine mincos und maxcos Funktionen? Unabhängig von der Eingabe einen Wert und nur diesen einen Wert zurückgeben. Das habe ich für mich etwas abgekürzt.
Kann man vielleicht etwas Zeit gewinnen, wenn man die Reihenfolge der Integrationen vertauscht?
Ich habe das Problem immer noch nicht zur Gänze verstanden (s.u.), aber warum sollte das helfen? Am Ende mußt Du immer noch addieren und da spielt es keine Rolle welches Teilintegral zuerst berechnet wird. Im Gegenteil: Auteilen mag es Dir einfacher machen zu testen, aber es kostet Zeit - wie kritisch das ist, verwag ich hier nicht zu sagen, wahrscheinlich nicht sehr kritisch.
Die Funktion ist sehr scharf.
Das mag sein. Wenn ich eine geschriebe Gleichung sehe, kann ich auch sofort entscheiden, ob ich erst potenzieren muß oder nicht. Bei Gleichungen im Code ist neige ich (und nicht nur ich allein) dazu das fehlerhaft umzusetzen - es sei denn ich setze "überflüssige" Klammern. Das mag für Mathematiker peinlich sein - ich stehe dazu.
Ich versuche mal alles einzubauen und insbesonder den code übersichtlicher zu gestallten.
Gut, achte besonders darauf alles mit bedeutungsschwangeren Kommentaren zu versehen, damit wir verstehen, wohin Du ganz genau willst ;-). (Ich bin auch kein Mathematiker.) Man kann nicht genug betonen, wie wichtig gut kommentierter Code ist - gerade, wenn man wie Du gezwungen ist, viele Abkürzungen zu verwenden, damit der Numerikteil nicht noch mehr aufgeblasen wird.

Noch eine Frage in diesem Zusammenhang: Woher komme die Zahlen in g()? Sie fallen etwas vom Himmel. Außerdem frage ich mich, warum sie so "kurz" sind: Viel ungenauer, als Deine Fehlertoleranz ;-).

Kennst Du f2py? Das ist Teil von numpy. Wenn es auf die Performance ankommt, kannst Du damit den Code Deines Chefs einbetten und von Python aus aufrufen, was Dir die ganze Power von Python bringt.

Gruß,
Christian
Heini
User
Beiträge: 21
Registriert: Donnerstag 12. April 2007, 10:48

Guten Morgen.

Also ich habe die neue Variante oben in das Fenster eingefügt damit hier nicht so viele lange Pogramme stehen die ohnehin nicht lesbar und auch nicht mehr aktuell sind.

- Die Doppelschleife verwende ich um die Intervalle zu unterteilen, da ich weiß dass sowohl tplquad als auch dblquad+gauß sehr schnell in sehr kleine Intervalle vorstoßen müssen.
In der jetztigen Variante unterteile ich das Intervall "manuell" in 2 Teile, da ich weiß dass sich die Schwierigkeiten für k und k4 nahe 0 abspielen. Die Integration bzgl. cos(theta) ist wirklich harmlos.
Das Ergebnis:
dblquad+gauss erechnet das Integral für i=0 recht hurtig. Für i=1 dauert es dann schon ziemlich lange und für i=2 schreibt er nur noch eine Warnung über das Ausreißen der "subdivision-grenze" von 50:
" Warning: The maximum number of subdivisions (50) has been achieved.
If increasing the limit yields no improvement it is advised to analyze
the integrand in order to determine the difficulties. If the position of a
local difficulty can be determined (singularity, discontinuity) one will
probably gain from splitting up the interval and calling the integrator
on the subranges. Perhaps a special-purpose integrator should be used."

tplquad liefert bei mir nach 15 Minuten auch für i=0 keinen Wert.

- Gibt es auch einen Unterschied zwischen 1/() und ()**-1

- tuple-unpacking? Was ist das und wie kann ich es benutzen? Bezieht sich das auf deine Beispieldefinition "foo"? Ich bin nicht sicher ob ich verstanden habe wie man das anwendet.

- Mit der Reihenfolge meine ich folgendes:
k in [0,100] und k4 in [-5,5]. Wenn er jetzt jedes mal für festes k4 die k-Integration ausführt und dann das betreffende Intervall erneut teilt und diese Integration dann nochmal macht, macht er dann am Ende mehr? Jedes mal wird er für die k-Integration sehr schnell beim Intervall [0,0.1] landen.
Ich weiß einfach nicht was dblquad und tplquad so wirklich tun.

- Mein Chef macht das wie folgt:
cos(theta)-Integration per Gauß. k, k4-Integration mit adaptive Simpson. Für letzteres unterteilt er die Intervalle in 0.1-Intervalle. Er benutzt dabei einen Fehler von 10**(-10). Damit bekommt er ein Ergebnis. Kannst du auch kruz ein Beispiel schreiben wie f2py anzuwenden ist?

- Ach ja, ich habe im ersten Versuch den Fehler auf 10**(-5) gesetzt, weil der Standardfehler bei 10**(-9) liegt und ich einfach wollte, dass er nicht so lange integriert.

- Die Zahlen in g() sind an numerisch berechnete Kurven gefittete Werte. Ich habe auch schon darüber nachgedacht, die tools zur Integration von samples zu benutzen, aber da machen mir wohl die unterschiedlichen Argumente in den g()-Funktionen zu schaffen, wenn sie in F() aufgerufen werden.

Gruß, Thomas
CM
User
Beiträge: 2464
Registriert: Sonntag 29. August 2004, 19:47
Kontaktdaten:

Hoi Thomas,
Heini hat geschrieben: Also ich habe die neue Variante oben in das Fenster eingefügt damit hier nicht so viele lange Pogramme stehen die ohnehin nicht lesbar und auch nicht mehr aktuell sind.
Ist schon ok, aber das hat auch den Nachteil, daß man ewig scrollen muß.
Heini hat geschrieben: - Die Doppelschleife verwende ich um die Intervalle zu unterteilen, da ich weiß dass sowohl tplquad als auch dblquad+gauß sehr schnell in sehr kleine Intervalle vorstoßen müssen.
Das ist natürlich ein guter Punkt.
Heini hat geschrieben: dblquad+gauss erechnet das Integral für i=0 recht hurtig. Für i=1 dauert es dann schon ziemlich lange und für i=2 schreibt er nur noch eine Warnung über das Ausreißen der "subdivision-grenze" von 50:
Sorry, dazu fällt mir nichts mehr ein. Ich mache das zwar ähnlich, aber ich benötige nicht diese Genauigkeit. Mein Vorschlag: Wende Dich ab hier an die scipy-user Mailingliste. Da ist die Chance groß, daß man Dir weiterhelfen kann. Versuche das Problem, inkl. Fragestellung, so genau als möglich zu beschreiben. Inzwischen sieht der Code durchaus so aus, daß die Leute nicht nur drüber herfallen, statt auf die Frage zu antworten ;-).
Heini hat geschrieben: - Gibt es auch einen Unterschied zwischen 1/() und ()**-1
So etwas kannst Du ganz einfach austesten:

Code: Alles auswählen

import time
def test():
    a = time.time()
    for x in xrange(1, 1000000):
        y = float(x)
        i = 1.0/y
    b = time.time()
    for x in xrange(1, 1000000):
        y = float(x)
        i = y**-1
    c = time.time()
    print b-a, c-b
    
if __name__ == "__main__":
    test()
(Billiger test, ich weiß, aber ich habe auch noch andere Dinge, die ich erledigen muß.)
Gibt bei mir: 1.72639012337 1.45548701286 Das ist (wenn auch mit anderen Werten), reproduzierbar. Wenn die Zahlen groß aussehen, so liegt das daran, daß auf meiner Maschine im Hintergrund Rechnungen laufen.
M. a. W.: Potenzieren ist etwas billiger.
Heini hat geschrieben: - tuple-unpacking? Was ist das und wie kann ich es benutzen? Bezieht sich das auf deine Beispieldefinition "foo"? Ich bin nicht sicher ob ich verstanden habe wie man das anwendet.
Siehe Turtorial.
(Meine Formulierung war nicht ganz hasenrein.)
Heini hat geschrieben: Ich weiß einfach nicht was dblquad und tplquad so wirklich tun.
Das Schöne an OpenSource: Man kann nachschauen. In Deinem Scipy-Download befindet sich: scipy/Lib/integrate/quadpack.h
Heini hat geschrieben: - Mein Chef macht das wie folgt:
cos(theta)-Integration per Gauß. k, k4-Integration mit adaptive Simpson. Für letzteres unterteilt er die Intervalle in 0.1-Intervalle. Er benutzt dabei einen Fehler von 10**(-10). Damit bekommt er ein Ergebnis. Kannst du auch kruz ein Beispiel schreiben wie f2py anzuwenden ist?
Klar: (Direkt aus der Anleitung ;-) .)

Code: Alles auswählen

C FILE: ARRAY.F
      SUBROUTINE FOO(A,N,M)
C
C     INCREMENT THE FIRST ROW AND DECREMENT THE FIRST COLUMN OF A
C
      INTEGER N,M,I,J
      REAL*8 A(N,M)
Cf2py intent(in,out,copy) a
Cf2py integer intent(hide),depend(a) :: n=shape(a,0), m=shape(a,1)
      DO J=1,M
         A(1,J) = A(1,J) + 1D0
      ENDDO
      DO I=1,N
         A(I,1) = A(I,1) - 1D0
      ENDDO
      END
C END OF FILE ARRAY.F
Am Besten Du ließt Dir die f2py-Anleitung durch - denn alles weitere führt an dieser Stelle zu weit.

Tut mir leid. Aber ich würde mich freuen, wenn Dir sonst jemand hier oder auf der scipy-Liste helfen kann. Wenn man Dir auf der scipy-Liste hilft, sei so gut und poste hier die Lösung. Am WE bin ich nicht da, aber ansonsten lese ich die scipy-Liste mit und kann ggf. auch noch helfend eingreifen, wenn es zu Mißverständnissen kommt - schließlich kenne ich inzwischen etwas besser der Problem.

Gruß,
Christian
Heini
User
Beiträge: 21
Registriert: Donnerstag 12. April 2007, 10:48

Also mit den zahlreichen Tips hast du mir schon sehr geholfen. Das Programm an sich ist ja schon besser geworden und für meinen allgemeinen Programmierstil war das wohl auch sehr wohltuend.
Ich werde dann mal versuchen die Mailing-Liste zu kontaktieren.

Gruß, Thomas
Goa
User
Beiträge: 18
Registriert: Dienstag 10. Oktober 2006, 11:00

Noch ein allgemeiner Tipp: Teilweise wiederholst du die gleichen Berechnungen tausende male. Ich denk da vor allem an die ganzen M**2, M**3...
ich würde das an deiner Stelle genau einmal direkt nach der Definition ausrechnen, das heisst du definierst dir neue konstanten a la M2, M3 (eventuell auch per dictionnary M[2],...) und ersetzt es in den Funktionen F und Z.
Den selben Trick kannst du natürlich auch innerhalb der Funktionen für k und k0 anwenden.
Keine Ahnung wieviel dir das an Rechenzeit erspart (würd mich allerdings interessieren, d.h. ein bischen Feedback wäre nett), aber kostet ja nix das zu implementieren.
CM
User
Beiträge: 2464
Registriert: Sonntag 29. August 2004, 19:47
Kontaktdaten:

Oh ja, das ist ein guter Punkt! Üblicherweise solltest so etwas im Kopf der Funktion machen. Den Lookup per Dictionary würde ich (aus dem Bauch heraus, ich habe keine Belege) nicht empfehlen: Der Performancehit ist nicht sehr groß, meistens kann man so etwas machen, aber hier ... ? Also bleibt, z.B.

Code: Alles auswählen

M2 = M**2
M3 = M**3
Na ja, das ist trivial, spart aber sicherlich Rechenzeit. Und weshalb ich das überhaupt nochmal schreibe: Ich frage mich, ob es nicht tatsächlich sinnvoll ist derartige Paramter global zu lassen: Sonst wird bei jedem Funktionsaufruf noch mehr mitgeschleppt, die Funktion schluckt einen Haufen Parameter, die auch nicht gerade übersichtlich aussehen werden und Errechnen der Potenzen im Kopf würde ebenfalls bedeuten diese Parameter stehts neu auszurechnen. Also: Einmal ausrechnen und dann global lassen (ausnahmsweise)? Was denkt ihr?

Gruß,
Christian
CM
User
Beiträge: 2464
Registriert: Sonntag 29. August 2004, 19:47
Kontaktdaten:

@Heini: Habe nichs auf der scipy-Liste gelesen? Ist das Problem gelöst? Habe ich die Diskussion völlig übersehen?

Gruß,
Christian
BlackJack

Wenn sich etwas nicht ändert, kann man eine globale Konstante daraus machen. Das ist an sich ja nicht unsauber, eher das neu binden von globalen Variablen.

Aus Geschwindigkeitsgründen könnte man diese Berechnungen oder Ergebnisse auch als Defaultargumente einführen. Die Namen sind dann schneller beim Zugriff weil sie nicht im Namensraum des Moduls ("dynamischer" Dictionary-Zugriff) sondern in der Funktion (Array-Zugriff mit festem Index) liegen.

Aber Grundsätzlich gilt immer noch das Motto: Erst sollte die Funktion das korrekte Ergebnis liefern können, dann kann man optimieren.
Antworten