Group Abstract Group Abstract

Message Boards Message Boards

1
|
7.2K Views
|
3 Replies
|
1 Total Like
View groups...
Share
Share this post:

How to deal with singularity in numerical integration?

Kat Z
Posted 10 years ago

I would like to find the numerical integration result for

 f=(\[Pi]^2 (-20000 Sqrt[x] (\[Pi]^2 + 40000000 x) + Sqrt[10] 
E^(\[Pi]^2/(40000000 x)) \[Pi]^(3/2) (\[Pi]^2 + 60000000 x) Erfc[\[Pi]/(2000 Sqrt[10] Sqrt[x])]) 
(-20000 Sqrt[-x + \[Gamma]] + Sqrt[10] E^(\[Pi]^2/(40000000 (-x + \[Gamma]))) \[Pi]^(3/2)
Erfc[\[Pi]/(2000 Sqrt[10]Sqrt[-x + \[Gamma]])]))/(320000000000000000000000 x^(7/2)Sqrt[-x + \[Gamma]])

When \[Gamma]=10^-5, I can do numerical integration from x=0 to \[Gamma]. However, I received an error message saying the integrand ... has evaluated to Overflow, Indeterminate or Infinity for all sampling points in the region with boundaries ..., when I tried to evaluate

NIntegrate[f3, {x, 0, \[Gamma]}, PrecisionGoal -> 10, WorkingPrecision -> 100]` for `\[Gamma]=10^-4

I have plotted the function f vs x and found that f goes to infinity when x is close to 0. I am wondering how to handle this singularity problem. Any help is much appreciated.
POSTED BY: Kat Z
3 Replies
Sort By:

Not sure what is the strange behavior of Plot and NIntegrate. But if you rescale your function to simplify:

 res[y_, g_] = f/(40000000 \[Pi]^2 ) /. {x -> y/40000000, \[Gamma] -> g/40000000} // FullSimplify;
 res[y, g] // TraditionalForm

enter image description here

and find 

y0[g_] = Limit[res[y, g], y -> 0] // FullSimplify

enter image description here

it looks like your function is not divergent at y = 0 in the whole range of g:

Plot[y0[g], {g, 0, 100}, PlotTheme -> "Detailed"]

enter image description here

I could be wrong though, did not have much time to dig in.
POSTED BY: Vitaliy Kaurov
Kat Z
Kat Z
Posted 10 years ago

Thanks Vitaliy for your help.

Why I said that there is a singular points at 0 is because I get the following figure when I plot the function f vs x with [Gamma]=10^-5:

\[Gamma] = 10^-5;
Plot[f, {x, 0, 10^-5}]

enter image description here

However, I found that when I plot the simplified equation you get with g=400, I get the following plot:

f1[y_]=res[y,400];
Plot[f1[y], {y, 0, 10^-5}]

enter image description here

It seems the function is actually oscillating. I am very confused now.

Any clue, please?
POSTED BY: Kat Z

I think this is due to some error / precision loss. You have a complicated function and need to first analyze it a bit analytically, maybe take it apart on simpler components to understand larger behavior. Then you could approach numerical integration. But nearest future I will not have time to play with this.
POSTED BY: Vitaliy Kaurov
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard