# Numerical Integration in Laplace domain

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. –  Ｊ. Ｍ. 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. –  Ｊ. Ｍ. 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

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... –  Ｊ. Ｍ. 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