Methode zur numerischen Integration eines schwierigen Schwingungsintegrals

25

Ich muss das folgende Integral numerisch auswerten:

0sinc(xr)rE(r)dr

wobei , und . Hier ist die modifizierte Bessel-Funktion der zweiten Art. In meinem speziellen Fall habe ich , und .xR+λ,κ,ν>0Kλ=0,00313κ=0,00825ν=0,33E(r)=r4(λκ2+r2)ν5/2Kν5/2(λκ2+r2)xR+λ,κ,ν>0Kλ=0.00313κ=0.00825ν=0.33

Ich benutze MATLAB und habe die eingebauten Funktionen integralund ausprobiert quadgk, wodurch ich viele Fehler erleide (siehe unten). Ich habe natürlich auch zahlreiche andere Dinge ausprobiert, wie das Integrieren nach Teilen und das Summieren von Integralen von zu .( k + 1 ) x πkxπ(k+1)xπ

Haben Sie Vorschläge, welche Methode ich als nächstes ausprobieren sollte?

UPDATE (hinzugefügte Fragen)
Ich habe den Artikel @Pedro gelesen, der mit diesem Artikel verlinkt ist, und ich glaube nicht, dass es zu schwer zu verstehen war. Ich habe jedoch ein paar Fragen:

  • Wäre es in Ordnung, als Basiselemente in der beschriebenen univariaten Levin-Methode zu verwenden?ψ kxkψk
  • Könnte ich stattdessen einfach eine Filon-Methode verwenden, da die Frequenz der Oszillationen fest ist?

Beispielcode
>> integral(@(r) sin(x*r).*sqrt(E(r)),0,Inf)
Warning: Reached the limit on the maximum number of intervals in use. Approximate
bound on error is 1.6e+07. The integral may not exist, or it may be difficult to
approximate numerically to the requested accuracy.
> In funfun\private\integralCalc>iterateScalarValued at 372
In funfun\private\integralCalc>vadapt at 133
In funfun\private\integralCalc at 84
In integral at 89

ans =

3.3197e+06

torbonde
quelle
Was ist in deinem Integral? x
Pedro
Jede positive, reelle Zahl. Ich habe gerade meinen Beitrag aktualisiert.
Torbonde
Wenn Sie Code und Fehler anzeigen könnten, ist es wahrscheinlich nicht zu schwierig, die meisten davon zu lösen. Bitte lesen Sie den Fehler zuerst sorgfältig durch und prüfen Sie, ob Sie ihn selbst verschwinden lassen können.
Dennis Jaheruddin
Ich werde heute einen Kommentar mit einigem Code und Fehlern abgeben. Oder morgen.
Torbonde
Okay, also habe ich vergessen. Aber jetzt habe ich meinen Beitrag mit einem Beispiel aktualisiert (ich habe das Integral durch explizite Berechnung von in zwei Teile geteilt ). sinc
Torbonde

Antworten:

12

Ich habe meinen eigenen Integrator geschrieben, quadccder mit Singularitäten wesentlich besser zurechtkommt als die Matlab-Integratoren und eine zuverlässigere Fehlerschätzung liefert.

Um es für Ihr Problem zu verwenden, habe ich Folgendes getan:

>> lambda = 0.00313; kappa = 0.00825; nu = 0.33;
>> x = 10;
>> E = @(r) r.^4.*(lambda*sqrt(kappa^2 + r.^2)).^(-nu-5/2) .* besselk(-nu-5/2,lambda*sqrt(kappa^2 + r.^2));
>> sincp = @(x) cos(x)./x - sin(x)./x.^2;
>> f = @(r) sincp(x*r) .* r .* sqrt( E(r) );

Die Funktion fist jetzt Ihr Integrand. Beachten Sie, dass ich gerade einen alten Wert zugewiesen habe x.

Um in eine unendliche Domäne zu integrieren, setze ich Variablen ein:

>> g = @(x) f ( tan ( pi / 2 * x ) ) .* ( 1 + tan ( pi * x / 2 ).^2 ) * pi / 2;

dh das Integrieren gvon 0 nach 1 sollte dasselbe sein wie das Integrieren fvon 0 nach . Verschiedene Transformationen können unterschiedliche Qualitätsergebnisse liefern: Mathematisch sollten alle Transformationen das gleiche Ergebnis liefern, aber verschiedene Transformationen können glattere oder leichter integrierbare s erzeugen .g

Ich rufe dann meinen eigenen Integrator an quadcc, der mit dem NaNs auf beiden Seiten umgehen kann :

>> [ int , err , npoints ] = quadcc( g , 0 , 1 , 1e-6 )
int =
  -1.9552e+06
err =
   1.6933e+07
npoints =
       20761

Beachten Sie, dass die Fehlerschätzung sehr groß ist, dh quadccnicht viel Vertrauen in die Ergebnisse hat. Bei Betrachtung der Funktion ist dies jedoch nicht verwunderlich, da sie bei Werten schwankt, die drei Größenordnungen über dem tatsächlichen Integral liegen. Die Verwendung einer anderen Intervalltransformation kann wiederum zu besseren Ergebnissen führen.

Möglicherweise möchten Sie auch spezifischere Methoden wie diese untersuchen . Es ist ein bisschen komplizierter, aber definitiv die richtige Methode für diese Art von Problem.

Pedro
quelle
Vielen Dank. Ich werde einen Blick auf die verschiedenen Methoden werfen. Für meine Zwecke muss der Fehler nicht so klein sein, wie es in Gleichung integral(1e-10, glaube ich) üblich ist, aber 1.7e + 07 ist immer noch sehr, sehr groß. Vielleicht hilft eine andere Transformation, wie Sie bereits erwähnt haben.
Torbonde
@ cimrg.joe: Beachten Sie, dass die Fehlerschätzung eine Schätzung des absoluten Fehlers ist, die unter anderem auf den maximalen absoluten Werten des Integranden basiert. In einigen extremen Fällen kann der zurückgegebene Wert tatsächlich ganz in Ordnung sein. Wenn Sie nach einer Genauigkeit von zehn Stellen suchen, empfehle ich nachdrücklich die Verwendung der Levin-Methode, die ich am Ende meines Beitrags erwähnt habe.
Pedro
Ich brauche vielleicht keine zehn Stellen Genauigkeit, aber ich glaube, ich brauche mindestens fünf. Kann deine Methode das produzieren?
Torbonde
Die Methode kann diese Genauigkeit für Ihr Integral nicht garantieren, da die Werte am rechten Ende des Intervalls mehrere Größenordnungen größer sind als das Integral selbst.
Pedro
11

Wie Pedro betont, sind Levin-Methoden die am besten etablierten Methoden für diese Art von Problemen.

Haben Sie Zugang zu Mathematica? Für dieses Problem erkennt und verwendet Mathematica sie standardmäßig:

In[1]:= e[r_] := 
 r^4 (l Sqrt[k^2 + r^2])^(-v - 5/2) BesselK[-v - 5/2, l Sqrt[k^2 + r^2]]

In[2]:= {l, k, v} = {0.00313, 0.00825, 0.33};

In[3]:= Block[{x = 10}, 
 NIntegrate[Sinc'[x r] r Sqrt[e[r]], {r, 0, \[Infinity]}, 
  PrecisionGoal -> 3]]

Out[3]= -112494.

Hier ist ein Diagramm über einen Wertebereich von x:

In[4]:= ListLinePlot[
 Table[NIntegrate[Sinc'[x r] r Sqrt[e[r]], {r, 0, \[Infinity]}, 
   PrecisionGoal -> 3], {x, .5, 10, 0.1}]]

Zeichnen Sie von x = 0,5 bis x = 10

Sie können die jeweilige anzuwendende Levin-Methode auch manuell angeben, was in diesem Fall zu einer leichten Leistungsverbesserung führen kann:

In[5]:= method = {"LevinRule", "Kernel" -> {Cos[r x], Sin[r x]}, 
   "DifferentialMatrix" -> {{0, -x}, {x, 0}}, 
   "Amplitude" -> {(
     3497.878840962873` Sqrt[(
      r^4 BesselK[-2.17`, 
        0.00313` Sqrt[
         0.00006806250000000001` + r^2]])/(0.00006806250000000001` + 
        r^2)^1.415`])/
     x, -((3497.878840962873` Sqrt[(
       r^4 BesselK[-2.17`, 
         0.00313` Sqrt[
          0.00006806250000000001` + r^2]])/(0.00006806250000000001` + 
         r^2)^1.415`])/(r x^2))}, "AdditiveTerm" -> 0};

In[6]:= Block[{x = 10}, 
 NIntegrate[Sinc'[x r] r Sqrt[e[r]], {r, 0, \[Infinity]}, 
  PrecisionGoal -> 3, Method -> method]]

Out[6]= -112495.

Weitere Informationen zu Levin-Methoden in Mathematica finden Sie in der Dokumentation .

Andrew Moylan
quelle
Leider habe ich keinen Zugang zu Mathematica - nur zu MATLAB. Ich werde meine Frage nur mit einigen zusätzlichen Fragen zu dem Artikel @Pedro aktualisieren, auf den verwiesen wird.
Torbonde
OK, wie Sie sagen, müssen Sie mit Matlab auskommen. Ich werde eine weitere Antwort hinzufügen.
Andrew Moylan
5

Wenn Sie keinen Zugang zu Mathematica haben, können Sie in Matlab eine Levin-Methode (oder eine andere spezialisierte Oszillationsmethode) schreiben, wie Pedro vorschlägt.

Verwenden Sie die Chebfun- Bibliothek für Matlab? Ich habe gerade erfahren, dass hier eine grundlegende Levin-Methode implementiert ist . Die Implementierung wurde von Olver (einem der Experten auf dem Gebiet der oszillierenden Quadratur) geschrieben. Es befasst sich nicht mit Singularitäten, adaptiver Unterteilung usw., aber es kann genau das sein, was Sie brauchen, um anzufangen.

Andrew Moylan
quelle
Ich habe überlegt, eine Levin-Methode selbst zu implementieren, bin mir aber noch nicht sicher, ob ich für die Herausforderung bereit bin. Ich denke, ich muss die Methode ein bisschen besser verstehen. Vielleicht könnte ich mit meinem Berater darüber sprechen. Der Grund, warum ich nach den Filon-Methoden gefragt habe, ist, dass sie einfacher zu implementieren sind. Und da ich extrem hohe Genauigkeit nicht brauchen, aber das ist Teil meiner Diplomarbeit, wiegt schwer in.
torbonde
Ich habe mir die Chebfun-Bibliothek (die beeindruckend ist) und das Levin-Integrationsbeispiel angesehen. Aber ich kann es nicht zum Laufen bringen. Ich habe hier eine Frage dazu gestellt .
Torbonde
0

Die von Pedro empfohlene Transformation ist eine großartige Idee. Haben Sie versucht, mit den Parametern in Matlabs "quadgk" -Funktion herumzuspielen? Mit der Pedro-Transformation können Sie beispielsweise Folgendes tun:
quadgk(f, 0.0+eps, 1.0-eps, 'AbsTol', eps, 'MaxIntervalCount', 100000)
Mit dieser Methode erhalte ich eine Lösung von:
-2184689.50220729
und es dauert nur 0,8 Sekunden (mit den oben genannten Werten: x = 10).
Walter Gander und Walter Gautschi haben eine Arbeit über adaptive Quadratur mit Matlab Code, den Sie auch verwenden können (Link hier )

Joel Vroom
quelle