0
$\begingroup$

I'm new to Mathematica. I have written a code for performing a convolution integral (as follows) but it seems to be giving out error messages:

My code is:

a[x_?NumericQ] := PDF[NormalDistribution[40, 2], x]
b[k_?NumericQ, x_?NumericQ] := 0.0026*Sin[1.27*k/x]^2
c[k_?NumericQ, x_?NumericQ] := {a[x]*b[k, x]}
d[k_?NumericQ] := NIntegrate[c[k, x], {x, 0, Infinity}]
Plot [d[k], {k, 0, 350}]

It gives the following error message multiple times:

Integrand c[0.00715,x] is not numerical at {x} = {124.67}

And the plot dosen't show up..! Does anyone have any suggestions on how to fix it ?

$\endgroup$
  • 1
    $\begingroup$ You can remove curly brackets from the definition of cto make it numeric instead of a list. $\endgroup$ – Shadowray May 20 '17 at 6:48
  • $\begingroup$ It takes a while, but it's doable: Plot[0.0026 NIntegrate[Sin[1.27 k/x]^2 Exp[-(x - 40)^2/8], {x, 0, ∞}]/(2 Sqrt[2 π]), {k, 0, 350}] $\endgroup$ – J. M. will be back soon May 20 '17 at 6:52
  • $\begingroup$ If you use Method -> "GaussKronrodRule" it's about ten times faster. $\endgroup$ – Michael E2 May 20 '17 at 13:04
  • $\begingroup$ Hmm, a related bug: Integrate[(0.0026 E^(-(1/8) (-40 + x)^2) Sin[(1.27 k)/x]^2)/(2 Sqrt[2 Pi]), {x, 0, Infinity}] returns a constant 0.0013 (independent of k) which does not agree with NIntegrate. It seems to be the limit as k approaches infinity. $\endgroup$ – Michael E2 May 20 '17 at 19:22
  • $\begingroup$ Thank you all...! @Shadowray, This process works for the Normal Distribution function, but yields an error for the function a[x_] := 0.00003*(x^2) (1 - x/79.2).....this time it's an NIntegrate error..saying that the integral fails to converge. Changing the limit of integration from infinity to 90 helps, but is there a way to do the infinity integral? $\endgroup$ – Epari shalini May 21 '17 at 4:52
2
$\begingroup$

The only function definition that needs to have its argument restricted to numeric values is that for d since that is the only function that uses numeric techniques.

Clear[a, b, c, d]

a[x_] = PDF[NormalDistribution[40, 2], x];

b[k_, x_] = 0.0026*Sin[1.27*k/x]^2 // Rationalize;

As @Shadowray pointed out in the comments, the output of c should be a scalar rather than a list.

c[k_, x_] = a[x]*b[k, x]

(*  (13*Sin[(127*k)/(100*x)]^2)/
   (E^((1/8)*(-40 + x)^2)*
      (10000*Sqrt[2*Pi]))  *)

Using suggestion by @MichaelE2 to use Method -> "GaussKronrodRule"

d[k_?NumericQ] := NIntegrate[c[k, x], {x, 0, Infinity},
  Method -> "GaussKronrodRule"]

Plot[d[k], {k, 0, 350}]

enter image description here

$\endgroup$
  • $\begingroup$ Thank you for your inputs! I tried to follow the same procedure for a different function: a[x_] := 0.00003*(x^2) (1 - x/79.2). The NIntegrate shows error: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. "NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in x near {x} = {8.16907*10^224}. NIntegrate obtained -1.815272304467489*10^55882 and 1.815272304467489`15.954589770191005*^55882 for the integral and error estimates" $\endgroup$ – Epari shalini May 21 '17 at 4:48
  • 2
    $\begingroup$ @Eparishalini - for the a given in your comment, the integral does not converge on {0, Infinity}. $\endgroup$ – Bob Hanlon May 21 '17 at 5:09

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.