SingularPointIntegrationRadius < Mathematica < Mathe-Software < Mathe < Vorhilfe
|
| Status: |
(Frage) beantwortet  | | Datum: | 12:50 Do 27.08.2009 | | Autor: | PaRu |
Ich habe ein Integral, dass ich nummerisch lösen möchte. Da die zu integrierende Funktion eine Polstelle enthält, möchte ich das ganze nach dem Principal Value (http://en.wikipedia.org/wiki/Cauchy_principal_value) lösen. Das ganze funktioniert auch prima, wenn ich Integrate[] benutzt, aber bei NIntegrate[] bekomme ich die Fehlermeldung, dass ich einen singularity point angeben muss. obwohl ich dies mache, bekomme ich trotzdem die Fehlermeldung. Kann mir jemand dabei weiterhelfen?
Anbei ein Minimalbeispiel...
Dateianhänge:
Anhang Nr. 1 (Typ: nb) [nicht öffentlich]
|
|
| |
|
Also es war ziemlich kompliziert ohne wirkliches Vorwissen reinzugehen, aber ich glaube, ich hab dein Problem gefunden:
Der Befehl, der dich zu deinem gewünschten Ziel bringt sieht folgender Maßen aus: (omega schreibe ich hier mal als M, damit keine Verwechslung entsteht)
| 1: | kk2[w_] :=
| | 2: | NIntegrate[(w f[x])/(x^2 - w^2), {x, -10, -w, w, 10},
| | 3: | Method -> {"PrincipalValue",
| | 4: | "SingularPointIntegrationRadius" -> w}];
| | 5: | kk3[w_] :=
| | 6: | NIntegrate[(w f[x])/(x^2 - w^2), {x, -10, w, -w, 10},
| | 7: | Method -> {"PrincipalValue",
| | 8: | "SingularPointIntegrationRadius" -> -w}];
| | 9: | p1 = Plot[kk2[M], {M, 0, 6}, Frame -> True,
| | 10: | GridLines -> Automatic];
| | 11: | p2 = Plot[kk3[M], {M, -6, 0}, Frame -> True,
| | 12: | GridLines -> Automatic];
| | 13: | Show[p1, p2, PlotRange -> {{-6, 6}, {-1.5, 1.5}}] |
So nun zur Erklärung:
Der Befehl NIntegrate erfordert, dass die Singularitäten (hier also [mm] $\omega$ [/mm] und [mm] $-\omega$) [/mm] in den Integrationsbereich mit eingetragen werden. Hier nun das erste Problem: Angenommen, [mm] $\omega$ [/mm] ist negativ, dann wäre die Reihenfolge nicht mehr gegeben, also muss man das ganze in den negativen und den positiven Bereich aufteilen.
Weiterhin muss der Radius angegeben werden. Der Ausdruck x^2 = M^2 liefert hier aber nur ein true/false zurück und keine Zahl! Der Radius muss direkt nach dem Pfeil als Wert angegeben werden.
Schreibt man nur [mm] $\omega$ [/mm] führt dies in negativen Bereichen zu unschönen Ergebnissen. Deswegen hier einmal mit [mm] $\omega$ [/mm] und [mm] $-\omega$. [/mm] Anschließend schnell die beiden Plots für den positiven und den negativen Bereich berechnet und dann in einem Plot zusammengefügt.
Bei mir sah das Ganze dann so aus:
[Dateianhang nicht öffentlich]
lg Sunny
Dateianhänge:
Anhang Nr. 1 (Typ: gif) [nicht öffentlich]
|
|
|
| |
|
| Status: |
(Frage) beantwortet | | Datum: | 10:20 So 30.08.2009 | | Autor: | PaRu |
danke ersteinaml für die schnelle antwort.
so wie es aussieht, scheint nun das ganze nicht auf grund der singularitypointradius-deklaration zu funktionieren, sondern wegen der exclusion von den polstellen (siehe anhang). es wundert mich, dass die definition des singularitypointradius keine auswirkung hat?
ich muss da nochmal nachhacken, da in meiner eigentlichen untersuchung mal wieder probleme beim lösen des integrals auftreten...
Dateianhänge:
Anhang Nr. 1 (Typ: nb) [nicht öffentlich]
|
|
|
| |
|
Na, ich ahnte doch, dass es eine einfache Lösung geben muss.
Es gibt ein eigentlich veraltetes Päckchen unter library.wolfram.com/infocenter/MathSource/6778 herunterzuladen und dann funktioniert:
| 1: | << CauchyPrincipalValue'
| | 2: | kk2a[w_] := CauchyPrincipalValue[w*(f[x]/(x^2 - w^2)),
| | 3: | {x, -10, {-Abs[w], Abs[w]/2}, {Abs[w], Abs[w]/2}, 10}]
| | 4: | Plot[kk2a[w], {w, -6, 6}] |
reibungslos.
[Dateianhang nicht öffentlich]
Noch'n Gruß,
Peter
P.S.: ich habe die Datei aus dem Netz in das Verzeichnis D:\Programme\Wolfram Research\Mathematica\7.0\AddOns\Applications gespeichert; dann klappt's auch so, wie oben angegeben (möglicherweise hast Du eine andere Partition, als "D:" - Hauptsache ist, dass Du die Datei in das in Deiner Konfiguration entsprechende Verzeichnis kopierst). So, jetzt habe ich aber fertig 
Dateianhänge:
Anhang Nr. 1 (Typ: jpg) [nicht öffentlich]
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig  | | Datum: | 11:15 So 30.08.2009 | | Autor: | Peter_Pein |
Hallöle PaRu,
100%ig sicher bin ich mir auch nicht, aber zunächst scheint es mir sinnlos zu sein, einer Option, die ..Radius heißt, Punkte zuzuordnen.
Ein nicht zufriedenstellendes Ergebnis bekomme ich mit | 1: |
| | 2: | kk2[w_] := NIntegrate[(w*f[x])/(x^2 - w^2),
| | 3: | Evaluate[Flatten[{x, -10, Through[{Min, Max}[{-w, w}]],10}]],
| | 4: | Method -> {"PrincipalValue",
| | 5: | "SingularPointIntegrationRadius" -> 1/10},
| | 6: | MaxRecursion -> 20];
| | 7: | Plot[kk2[\[Omega]], {\[Omega], -6, 6}, Frame -> True, GridLines -> Automatic] |
Irgendwie stehe ich heute auf'm Schlauch und komme nicht weiter.
Aber ein Tipp zum Plotten einer exakt integrierbaren Funktion:
| 1: | Plot[kk1[\[Omega]], {\[Omega], -6, 6}, Frame -> True,
| | 2: | GridLines -> Automatic] // Timing |
ergibt auf meiner Mühle 105 und ein paar zerdrückte Sekunden. Erst | 1: | tmp = Assuming[Im[\[Omega]] == 0, FullSimplify[PiecewiseExpand[
| | 2: | kk1[\[Omega]]]]] |
und dann tmp plotten ergibt 0.02 Sekunden.
Aber ich werde noch versuchen dem Hauptwert in der numerischen Variante auf die Schliche zu kommen.
Gruß,
Peter
P.S.: Meine Funktionentheorievorlesungen sind zwar schon über zwanzig Jahre her, aber soweit ich mich erinnere, ist es bei der Integration um Singularitäten herum wichtig, dass der Integrationsweg nur eine Singularität umschließt und es ist egal, ob der Radius darum nun ein Millionstel ist oder 5 Milliarden (um mal konkrete Zahlen zu nennen).
Die Cauchy'sche Integralformel ist Dir ein Begriff, nehme ich an.
|
|
|
| |
|
| Status: |
(Frage) überfällig  | | Datum: | 11:43 Mo 31.08.2009 | | Autor: | PaRu |
Hallo Peter,
danke für die tipps.
die veraltete library hatte ich auch gefunden, doch leider kann ich sie nicht laden, obwohl ich sie sowohl im dem von dir genannten verzeichnis und dem verzeichnis mit dem quellcode hinterlegt hatte? (siehe anhang; ich benutze übrigens auch mathematica 7.0)
der vorschlag mit Assuming[] bezüglich der plot-rechenzeit hat leider bei mir keine verbesserung gebracht.
so, und nun nochmal zum eigentlichen problem: im anhang kann man sehen, dass bei meiner eigentlichen untersuchung, ein falsches ergebnis entsteht, wenn man nach der bisherigen variante integriert und plotet. wenn ich aber das integral wie bisher löse und mittels einer liste und ListPlot[] den graphen darstellen lasse, erhalte ich ein brauchbares ergebnis???
Gruß Patrick
ps: eigentlich hast du mit dem integrationsradius recht, aber evtl. spielt dieser bei der nummerischen integration doch eine rolle. wer weiss, welchen algorithmus mathematica zum lösen benutzt...
Dateianhänge:
Anhang Nr. 1 (Typ: nb) [nicht öffentlich]
|
|
|
| |
|
| Status: |
(Mitteilung) Reaktion unnötig | | Datum: | 12:20 Do 01.10.2009 | | Autor: | matux |
$MATUXTEXT(ueberfaellige_frage)
|
|
|
|
