# 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. Commented Jun 6, 2012 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... Commented Jun 6, 2012 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. Commented Jun 6, 2012 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. Commented Jun 6, 2012 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? Commented Sep 10, 2012 at 0:53

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? Commented Jun 7, 2012 at 13:28
• @IstvánZachar, it's correct. Nothing changed. Commented Jun 7, 2012 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… Commented Jun 8, 2012 at 13:35
• @IstvánZachar that gave me 15 points less ... Commented Sep 10, 2012 at 0:50

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... Commented Jun 6, 2012 at 15:45
• @Bruno: then you should try the suggestion I gave in my previous comment... Commented Jun 6, 2012 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… Commented Jun 6, 2012 at 17:55