Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
    
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
add comment

2 Answers

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.

share|improve this answer
    
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
add comment

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

]
share|improve this answer
    
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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