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?

Improve this question
$\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 missing motivation ♦
    Commented Jun 6, 2012 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$
    – Bruno Rangel
    Commented Jun 6, 2012 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 missing motivation ♦
    Commented Jun 6, 2012 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$
    – Peter Breitfeld
    Commented Jun 6, 2012 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$
    – Rolf Mertig
    Commented Sep 10, 2012 at 0:53

2 Answers 2

Reset to default
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}]

]
Improve this answer
$\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$
    – István Zachar
    Commented Jun 7, 2012 at 13:28
  • $\begingroup$ @IstvánZachar, it's correct. Nothing changed. $\endgroup$
    – Bruno Rangel
    Commented Jun 7, 2012 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$
    – Bruno Rangel
    Commented Jun 8, 2012 at 13:35
  • $\begingroup$ @IstvánZachar that gave me 15 points less ... $\endgroup$
    – Rolf Mertig
    Commented Sep 10, 2012 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.

Improve this answer
$\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$
    – Bruno Rangel
    Commented Jun 6, 2012 at 15:45
  • $\begingroup$ @Bruno: then you should try the suggestion I gave in my previous comment... $\endgroup$
    – J. M.'s missing motivation ♦
    Commented Jun 6, 2012 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$
    – Bruno Rangel
    Commented Jun 6, 2012 at 17:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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