Skip to main content

plotting - How to plot the result of this singular integral?


Please I open a new post here after this one : https://mathematica.stackexchange.com/a/59203/10158



Now I want to plot the function $f(a,b)$ as a function of $b$ for different values of $a$ : $a=0.5$ , $a=10$ and $a=100$


$f(a,b)=1-\dfrac{a \, b^2}{K_{2}(a)}\dfrac{\partial^2}{\partial a^2}\Bigg[\dfrac{1}{a}\dfrac{1}{\sqrt{1-b^2}}\displaystyle{\int_1^{\infty }\dfrac{e^{-ax}}{\sqrt{x^2-1} (1+x\sqrt{1-b^2}) (x-1/\sqrt{1-b^2})} \, dx}\Bigg]$


where $K_{2}(a)$ is the modified Bessel function of the second kind and $\dfrac{\partial^2} {\partial a^2}$ is the $2^{nd}$ derivative with respect to a.


Mathematica input form:


1 - ((a b^2/BesselK[2, a]) D[(1/a) (1/Sqrt[1 - b^2]) Integrate[E^-ax/(Sqrt[x^2 - 1] (1 + x Sqrt[1 - b^2]) (x - 1/Sqrt[1 - b^2])), {x,1, Infinity}], {a, 2}])

how do it please?


Thank's.



Answer



First I assume that the OP wants to compute the Cauchy principal value of the integral, since the integrand has a simple pole in the interior of the interval of integration. Dealing with such singularities is described in the tutorial NIntegrate Integration Strategies.



I let the derivative be computed in the process of defining f. The symbols a and b are highlighted in red by the syntax highlighter because they appear as function arguments and as symbols being blocked. It indicates a warning of a potential scoping conflict. But Block is evaluated before the definition of f is set, so the warning may be ignored in this case. The differentiation results in an expression that contains a linear combination of three NIntegrate integrals. I combine these, more or less manually, with pattern replacements that implement the linearity properties of the integral. This reduces computation time by about 40%.


ClearAll[f];

f[a_, b_] /; a > 0 && 0 < b < 1 := 
 Evaluate[Block[{a, b, NIntegrate}, 
   1 - ((a b^2/BesselK[2, a]) *
       D[(1/a) (1/Sqrt[1 - b^2]) *
          NIntegrate[Exp[-a x]/(Sqrt[x^2 - 1] (1 + x Sqrt[1 - b^2]) (x - 1/Sqrt[1 - b^2])),
           {x, 1, 1/Sqrt[1 - b^2], Infinity}, 
           Method -> "PrincipalValue"],

        {a, 2}]) /. 
    HoldPattern[Times[NIntegrate[int_, args__], coeff__]] :>
      NIntegrate[Times[int, coeff], args] /. 
    t : Plus[_NIntegrate, __NIntegrate] :> 
      NIntegrate[First /@ t, Sequence @@ t[[1, 2 ;;]]]
   ]]

For the plot, I set a few options for NIntegrate suitable for plotting. Somehow some of the coefficients get set to MachinePrecision which causes a precision warning; therefore I turned the warning off.


Module[{opts = Options[NIntegrate], plot, offQ},
 SetOptions[NIntegrate, {WorkingPrecision -> 12, AccuracyGoal -> 3, MaxRecursion -> 12}];

 offQ = Head[NIntegrate::precw] === $Off;
 Off[NIntegrate::precw];

 plot = Plot[
   Evaluate@Table[f[a, b], {a, {0.5`20, 10, 100}}], {b, 0, 1}, 
   PlotPoints -> 15, WorkingPrecision -> 12];

 SetOptions[NIntegrate, opts];
 If[! offQ, On[NIntegrate::precw]];
 plot

 ]

Mathematica graphics


Evaluation is still rather slow, the plot above taking about 28-29 seconds.


Comments

Post a Comment

Popular posts from this blog

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...
Post a Comment
Read more

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...
Post a Comment
Read more

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more