I need to calculate two different integrals containing a Bessel function in the Laplace domain. I have tried different kinds of quadrature but didn’t have any luck. I don’t know how to help Mathematica along to treat the Laplace variable (s) in the equations:

@ RolfMertig, thanks to the details. I have successfully added the inversion algorithm for time domain. But the results do not look very good and take a very long time to compute… This is really bugging me out! And I still couldn’t find a way to compute this efficiently…

eq1[(s_)?NumericQ] := NIntegrate[BesselK[0, Sqrt[u^2 + s]], {u, 0, 40}];

eq2[(s_)?NumericQ] := NIntegrate[BesselK[0, Sqrt[u^2 + s]], {u, 0, -10}];

eq3[(s_)?NumericQ, (n_)?NumericQ] :=

NIntegrate[BesselK[0, Sqrt[(-Log[t])^2 + ((n^2*Pi^2)/100 + s)]]*

Cos[n*Pi*((50 + (15 – Log[t]/Sqrt[(n^2*Pi^2)/100 + s])*

Cot[30*Degree])/100)]*(1/t), {t, 0, -10*Sqrt[(n^2*Pi^2)/100 + s]}];

eq4[(s_)?NumericQ, (n_)?NumericQ] :=

NIntegrate[BesselK[0, Sqrt[(-Log[t])^2 + ((n^2*Pi^2)/100 + s)]]*

Cos[n*Pi*((50 + (15 – Log[t]/Sqrt[(n^2*Pi^2)/100 + s])*

Cot[30*Degree])/100)]*(1/t), {t, 0, 40*Sqrt[(n^2*Pi^2)/100 + s]}];

PD[s_] = (1/(s^(3/2)*100*Sin[30*Degree]))*(eq1[s] – eq2[s]) +

(2/(s*100*Sin[30*Degree]))*Sum[(1/Sqrt[(n^2*Pi^2)/100 + s])*

Cos[n*(Pi/2)]*(eq3[s, n] – eq4[s, n]), {n, 1, 10}];

SetOptions[NIntegrate, Method -> “ClenshawCurtisRule”,

WorkingPrecision -> 5, MaxRecursion -> 12, AccuracyGoal -> 2];

Ni = 8;

V[i_, NN_] = (-1)^(i + NN/2)*Sum[(k^(1 + NN/2)*(2*k)!)/((i – k)!*

k!^2*(-i + 2*k)!*(-k + NN/2)!), {k, Floor[(1 + i)/2], Min[i, NN/2]}];

PD1[tD_] = (Log[2]*Sum[PD[(i*Log[2])/tD]*V[i, Ni], {i, 1, Ni}])/tD;

LogLogPlot[{PD1[y], y*Derivative[1][PD1][y]}, {y, 1, 100},

PlotStyle -> {{Black}, {Dashed, Black}}, Frame -> True]

I’m not sure if I’m missing something obvious or computing these equations in a completely rough way. Can anyone provide an insight?

=================

Is there any reason why you absolutely have to use Clenshaw-Curtis? Your eq1 seems to be fine, but I’d have done something like eq1[s_?NumericQ] := NIntegrate[BesselK[0, Sqrt[u^2 + s]], {u, 0, -10}, Method -> “ClenshawCurtisRule”] myself.

– J. M.♦

Jun 6 ’12 at 15:29

@J.M., the reason for Clenshaw-Curtis is that I found a reference that suggests the chebyshev polynomials as integration procedure for these computations. But it seems that no other method can work these equations as well…

– Bruno Rangel

Jun 6 ’12 at 15:38

You should have edited your question instead of posting a new answer. In any event: Your PD[] has calls to eq1 and others when they should be eq1[s]; also your functions being numeric preclude constructions like PD1′[y]; you will have to implement a separate method for the derivative if you need it.

– J. M.♦

Jun 6 ’12 at 18:17

eq3and eq4 have “n” it it, so you should use eq3[s_?NumericQ,n_?NumericQ] for these and call them in PD with 2 parameters. And you should use PD[s_]:=…. Than you can call e.g. PD[2] and get a value (0.0182218 + 0.0011354 I). But the integrals seem to convert badly, because you’ll get a bunch of messages, so you can’t trust this value.

– Peter Breitfeld

Jun 6 ’12 at 21:57

@BrunoRangel As J.M. already pointed out: your Derivative[1][PD1][y] does not make sense. You need to program it (and the derivatives for eq1, eq2, eq3, eq4). Maybe you can fix this?

– Rolf Mertig

Sep 10 ’12 at 0:53

=================

2 Answers

2

=================

You have 2 variables in your equations, but you integrate over one only. So the integrand will never be numeric (always a “s” in eq1). This is why NIntegrate yields the error. It you will give an numerical value for s, e.g, s=2 then NIntegrate will give you a numerical result.

NIntegrate always complains, if the result is not a number.

Thanks Peter, but the problem is that I don’t think that I can give the Laplace variable “s” a value. Because eventually I’ll need it for a numerical inversion to time domain…

– Bruno Rangel

Jun 6 ’12 at 15:45

@Bruno: then you should try the suggestion I gave in my previous comment…

– J. M.♦

Jun 6 ’12 at 15:51

Hey J.M.: Thanks for the suggestion! I have used, but Mathematica still can’t give me an expression to put into the inversion algorithm… I’m posting my codes bellow, so you can have a clear idea of the problem. I still can’t figure it out what it’s going wrong…

– Bruno Rangel

Jun 6 ’12 at 17:55

This will calculate hundred PD in a few minutes, depending on your computer:

AbsoluteTiming[

eq1[(s_)?NumericQ] := (Print[“eq1[“,s//InputForm,”]”]; NIntegrate[BesselK[0, Sqrt[u^2 + s]], {u, 0, 40}]);

eq2[(s_)?NumericQ] := (Print[“eq2[“,s//InputForm,”]”]; NIntegrate[BesselK[0, Sqrt[u^2 + s]], {u, 0, -10}]);

eq3[(s_)?NumericQ, (n_)?NumericQ] := (Print[“eq3[“,s//InputForm,”,”,n, “]”];

NIntegrate[BesselK[0, Sqrt[(-Log[t])^2 + ((n^2*Pi^2)/100 + s)]]*

Cos[n*Pi*((50 + (15 – Log[t]/Sqrt[(n^2*Pi^2)/100 + s])*

Cot[30*Degree])/100)]*(1/t), {t, 0, -10*Sqrt[(n^2*Pi^2)/100 + s]}]

);

eq4[(s_)?NumericQ, (n_)?NumericQ] := (Print[“eq4[“,s//InputForm,”,”,n,”]”];

NIntegrate[BesselK[0, Sqrt[(-Log[t])^2 + ((n^2*Pi^2)/100 + s)]]*

Cos[n*Pi*((50 + (15 – Log[t]/Sqrt[(n^2*Pi^2)/100 + s])*

Cot[30*Degree])/100)]*(1/t), {t, 0, 40*Sqrt[(n^2*Pi^2)/100 + s]}]

);

PD[s_] := (1/(s^(3/2)*100*Sin[30*Degree]))*(eq1[s] – eq2[s]) +

(2/(s*100*Sin[30*Degree]))*Sum[(1/Sqrt[(n^2*Pi^2)/100 + s])*

Cos[n*(Pi/2)]*(eq3[s, n] – eq4[s, n]), {n, 1, 10}];

ParallelEvaluate@

SetOptions[NIntegrate, Method -> “ClenshawCurtisRule”, WorkingPrecision -> 80, MaxRecursion -> 30, AccuracyGoal -> 2];

ParallelTable[PD[i],{i,100}]

]

Sorry, I might have mixed up this answer by accepting a suggested review by @Bruno Rangel. How the review improved the answer was unclear so I rolled it back. Could you please make sure that no harm was done?

– István Zachar

Jun 7 ’12 at 13:28

@IstvánZachar, it’s correct. Nothing changed.

– Bruno Rangel

Jun 7 ’12 at 13:43

@RolfMertig, thanks to the details. I have successfully added the inversion algorithm for time domain. But the results do not look very good and take a very long time to compute… This is really bugging me out! And I still couldn’t find a way to compute this efficiently…

– Bruno Rangel

Jun 8 ’12 at 13:35

@IstvánZachar that gave me 15 points less …

– Rolf Mertig

Sep 10 ’12 at 0:50