1
$\begingroup$

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?

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – J. M.'s torpor
    Jun 6 '12 at 15:29
  • $\begingroup$ @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... $\endgroup$ Jun 6 '12 at 15:38
  • $\begingroup$ 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. $\endgroup$
    – J. M.'s torpor
    Jun 6 '12 at 18:17
  • $\begingroup$ 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. $\endgroup$ Jun 6 '12 at 21:57
  • $\begingroup$ @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? $\endgroup$ Sep 10 '12 at 0:53
4
$\begingroup$

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}]

]
$\endgroup$
4
  • $\begingroup$ 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? $\endgroup$ Jun 7 '12 at 13:28
  • $\begingroup$ @IstvánZachar, it's correct. Nothing changed. $\endgroup$ Jun 7 '12 at 13:43
  • $\begingroup$ @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… $\endgroup$ Jun 8 '12 at 13:35
  • $\begingroup$ @IstvánZachar that gave me 15 points less ... $\endgroup$ Sep 10 '12 at 0:50
3
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ 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... $\endgroup$ Jun 6 '12 at 15:45
  • $\begingroup$ @Bruno: then you should try the suggestion I gave in my previous comment... $\endgroup$
    – J. M.'s torpor
    Jun 6 '12 at 15:51
  • $\begingroup$ 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… $\endgroup$ Jun 6 '12 at 17:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.