Skip to main content

numerical integration - Strategies to solve an oscillatory integrand only known numerically


I have an integrand that looks like this:


enter image description here


the details of computation are complicated but I only know the integrand numerically (I use NDSolve to solve second order ODE). The integrand is not simply the solution of my ODE either; calling the two solutions of my ODE osc1[s], osc2[s] then schematically the integrand I have looks something like exp(-is)[g(s)osc1[s]osc2*[C-s]+f(s)osc2[s]osc2*[C-s]]. The exp bit is only very slowly oscillating over my integration range, it is really osc1,osc2 that give wild oscillation, as a certain parameter they depend on gets larger.


More explicitely


rstar[r_] := r + 2 M Log[r/(2 M) - 1];
M=1;
rinf=10000;
rH = 200001/100000;
r0 = 10;

wp=40;
ac=wp-8;
\[Lambda][l_] = l (l + 1);

eq[\[Omega]_,l_] := \[CapitalPhi]''[r] + (2 (r - M))/(
r (r - 2 M)) \[CapitalPhi]'[
r] + ((\[Omega]^2 r^2)/(r - 2 M)^2 - \[Lambda][l]/(
r (r - 2 M))) \[CapitalPhi][r] == 0;
init=-0.0000894423075560122420468703835499 + 
0.0000447222944185058822813688948339 I;

dinit=-4.464175354293244250869336196691640386266791`30.*^-6 - 
8.950483248390306670770345406047835993931665`30.*^-6 I;

osc1 := \[CapitalPhi] /. 
Block[{$MaxExtraPrecision = 100}, 
NDSolve[{eq[1/10, 1], \[CapitalPhi][rinf] == 
init, \[CapitalPhi]'[rinf] == dinit}, \[CapitalPhi], {r, r0, 
rinf}, WorkingPrecision -> wp, AccuracyGoal -> ac, 
MaxSteps -> \[Infinity]]][[1]];


osc2 is obtained simiarly. Note these are for non problematic params and it will run quick quickly, and not be too badly behaved.


The problem I have is that I only know the integrand to maybe 6-12 digits of precision (dp), depending on the parameters. This is computing the NDSolve with a WorkingPrecision of 50-60, AccuracyGoal->42-52 and it takes around 2 hrs. I want to integrate this with NIntegrate, but when my parameters are large (and the oscillation is very high) I usually only know the integrand around the 6 dp end of scale, and NIntegrate wants a greater WorkingPrecision than this otherwise it complains (since oscillation is also getting very large).


I can force it to do the integral by making the WorkingPrecision higher, but I think this is cheating if I don't believe my integrand any higher than 6 dp?


The only ideas I've had so far are to try different rules. Are there any rules people would recommend for doing such oscillatory integrands? So far I've tried "LevinRule", "ClenshawCurtisRule", "GaussKronrodRule" but none seem to compute it any quicker than just the default. They all agree up to a reasonable number of dp, so no idea if I should just stick to the default, or if there is something better one could do with such an integrand. Speed is not a concern just accuracy.


UPDATE


Let's say I managed to split my integral into a few different integrals. First give the definitions:


vbar[tau_?
NumericQ] := (4 M) ((tau/tauh)^(1/3) + 1) Exp[-(tau/tauh)^(1/3) + 
 1/2 (tau/tauh)^(2/3) - 1/3 (tau/tauh)];
ubar[tau_?

NumericQ] := -(4 M) ((tau/tauh)^(1/3) - 1) Exp[(tau/tauh)^(1/3) + 
 1/2 (tau/tauh)^(2/3) + 1/3 (tau/tauh)];
rtau[tau_?NumericQ] := (2 M) (tau/tauh)^(2/3);

in addition to those made above, then I think I can give my integral as a sum of integrands that look like this


Exp[-I s] (ubar[tau_f - s])^(-i 4/10)Exp[+i 1/10 rstar[tau_f - s]]osc1[rtau[tauf-s]]*

here tau_f constant. The first part is an amplitude, the osc1 satisfies the linear ODE given above. I think this has Levin potential if I can work out how to input the LevinRules given the above second order ODE? (Here and in the above I fix my parameters the ODE depends on to (1/10,1) to simplify giving the ICs but I don't that detracts from the main problem). Would need to work out what the Kernal is from the ODE above.




Comments

Post a Comment

Popular posts from this blog

plotting - Filling between two spheres in SphericalPlot3D

Manipulate[ SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, Mesh -> None, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], {n, 0, 1}] I cant' seem to be able to make a filling between two spheres. I've already tried the obvious Filling -> {1 -> {2}} but Mathematica doesn't seem to like that option. Is there any easy way around this or ... Answer There is no built-in filling in SphericalPlot3D . One option is to use ParametricPlot3D to draw the surfaces between the two shells: Manipulate[ Show[SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], ParametricPlot3D[{ r {Sin[t] Cos[1.5 Pi], Sin[t] Sin[1.5 Pi], Cos[t]}, r {Sin[t] Cos[0 Pi], Sin[t] Sin[0 Pi], Cos[t]}}, {r, 1, 2 - n}, {t, 0, Pi}, PlotStyle -> Yellow, Mesh -> {2, 15}]], {n, 0, 1}]
Post a Comment
Read more

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]
Post a Comment
Read more

plotting - Mathematica: 3D plot based on combined 2D graphs

I have several sigmoidal fits to 3 different datasets, with mean fit predictions plus the 95% confidence limits (not symmetrical around the mean) and the actual data. I would now like to show these different 2D plots projected in 3D as in but then using proper perspective. In the link here they give some solutions to combine the plots using isometric perspective, but I would like to use proper 3 point perspective. Any thoughts? Also any way to show the mean points per time point for each series plus or minus the standard error on the mean would be cool too, either using points+vertical bars, or using spheres plus tubes. Below are some test data and the fit function I am using. Note that I am working on a logit(proportion) scale and that the final vertical scale is Log10(percentage). (* some test data *) data = Table[Null, {i, 4}]; data[[1]] = {{1, -5.8}, {2, -5.4}, {3, -0.8}, {4, -0.2}, {5, 4.6}, {1, -6.4}, {2, -5.6}, {3, -0.7}, {4, 0.04}, {5, 1.0}, {1, -6.8}, {2, -4.7}, {3, -1.
Post a Comment
Read more