Trying to crack a tough integral (Math nerds, help!)

[McC]

I've got this integral that I'm trying to solve.


G, O, F, a, and b are all unique (but changeable) constants. The only variable is r.

Ideally, I'd like to get it into such a form that I can compute P - P[a] directly.

I was up until around 3am trying to crack it this morning, with no luck. I just don't remember enough from my high school calculus classes.

[Simon_Jester]

I've got this integral that I'm trying to solve.


G, O, F, a, and b are all unique (but changeable) constants. The only variable is r.

Ideally, I'd like to get it into such a form that I can compute P - P[a] directly.

I was up until around 3am trying to crack it this morning, with no luck. I just don't remember enough from my high school calculus classes.

Wait... that's 0.5 to the power of the quadratic function?

I could swear I've seen this one in books of integral tables. Can you find a copy of the CRC Math handbook?

[Hawkwings]

Integrals.com does these with minimal effort.

http://integrals.wolfram.com/index.jsp

and the input I used is: 0.5^(((x-O)/f)^2)
Since G is just a multiplier, I left it out.

Code: Select all

                        2
             ((x - o)/f)
Integrate[0.5            , x] == 

                        (0. + 0.832555 I) (o - 1. x)
(0. + 1.06447 I) f Erfi[----------------------------]
                                     f

Only problem is, can you interpret that result?

[McC]

Only problem is, can you interpret that result?

No. This was actually the first place I went, but I am completely unfamiliar with "Erfi".

I also have a vague suspicion that it's only presenting that solution because it's trying to integrate from -inf to +inf.

The particular area I'm interested in is bounded by r >= O, so that the numerator in that fraction is either 0 or positive.

[Simon_Jester]

Yeah, that's a problem with computer programs for solving integrals; the solutions are not consistently helpful. It's why I prefer to use integral tables when I can.

Could you explain the underlying physical problem that requires you to do the integral? I'm curious.

[Xenophobe3691]

Wait. Just do substitution! 0.5 can go outside the integral (as it's a constant), and replace (r-O)/F with a simpler variable, integrate, then substitute back!

[McC]

Could you explain the underlying physical problem that requires you to do the integral? I'm curious.

It's just a curve on a graph that I want to find the area underneath.

Wait. Just do substitution! 0.5 can go outside the integral (as it's a constant), and replace (r-O)/F with a simpler variable, integrate, then substitute back!

0.5 is a constant, but it's being raised to a power, so it can't "just" go outside, so far as I know.

[Simon_Jester]

Could you explain the underlying physical problem that requires you to do the integral? I'm curious.

It's just a curve on a graph that I want to find the area underneath.

Yes, but WHY do you want the area under the curve? It matters. Do you need the area itself, for some specific set of known cases, or do you need the function for the area under the curve?

If you are trying to compute the area under a known curve, your best bet is to use an approximation method. There are a number of them. The best extremely simple one is the midpoint rule, and the best relatively simple one is Simpson's rule.

Ideally, you'd program that into a computer program such as Matlab or Mathematica. Such programs can calculate the are under a curve to many decimal places without having to actually do the integral. This is a very common way of avoiding nasty integrals- solve them computationally. Using computational methods to solve an integral by hand is far more time-consuming, but it is at least possible.

But that only works if you need the area under the curve but not the function for the area under the curve.
_______

If this is something like a homework problem, it's likely that you'll have to solve the problem analytically. Which means giving a function for the area under the curve. I strongly suggest that you look in libraries for reference books that contain integral tables, such as the CRC mathematics handbook.
_______

As for solving the integral by normal means, I have one idea. It's difficult, but just might work:

OK. Let's consider the following: you are integrating an expression of the form A^(((r-O)/F)^2) over r.

What happens if we define a variable x such that

x = ((r-O)/F)^2
________

In that case, there exists a function that describes r(x). It'll be messy, but hopefully you can at least prove that it has only one real root, depending on the value of O.

IF you can calculate r(x), you can get dr/dx. In which case you can say:

[Integral of]f(r)dr = [Integral of][f(x)(dr/dx)dx].

But dr/dx will be some (messy) polynomial or square root or something, while f(x) will simply be 0.5^x. The integral of 0.5^x is fairly straightforward, and hopefully you can use the product rule to hash out the mess you get when you include the polynomial.

I do not recommend this. But it's the only thing I can think of without actually spending several hours on the problem, which I am not inclined to do right now.

[Surlethe]

Here are a couple of steps to solve the integral.
Spoiler

First, change variables to p = (r - O)/F, and dr = Fdp, so that your integrand becomes F 0.5^[p^2] dp. The integral limits become α = (a - O)/F and β = (b - O)/F.

Spoiler

Next, change the base of the exponent to e: 0.5^[p^2] = exp[p^2 log(0.5)].

Spoiler

Now your integral has the form P = FG&integral;_{α}^{β} exp[p^2 log(0.5)] dp. I'm afraid it appears to have no simple antiderivative or closed form solution; the list of integrals on Wikipedia says the closed-form answer would be complex: FGsqrt[pi/(-4olog(0.5))]erf(sqrt[plog(0.5)]), where "erf" is the error function (look it up on Wikipedia).

Edit:Spoiler

Program a routine to do it numerically that takes your bounds and a fineness measure as inputs.

[Surlethe]

What happens if we define a variable x such that

x = ((r-O)/F)^2

I can tell you right away that won't work. Can you figure out how to quickly see why? Spoiler

The variable r only resides in the exponent, so if you replace r^2 with x you'll have dr = dx/(2r), but there is no 1/r term in the integrand.

In that case, there exists a function that describes r(x). It'll be messy, but hopefully you can at least prove that it has only one real root, depending on the value of O

I'm not sure why you want to find roots. The important thing here is to calculate dr and relate it to dx, which is not difficult to do implicitly: Spoiler

dx/dr = 2[(r-O)/F^2], so dr = F^2 dx / [2(r-O)]. Here you run into the problem to which I alluded above.

[McC]

Here are a couple of steps to solve the integral.
Spoiler

First, change variables to p = (r - O)/F, and dr = Fdp, so that your integrand becomes F 0.5^[p^2] dp. The integral limits become α = (a - O)/F and β = (b - O)/F.

Spoiler

Next, change the base of the exponent to e: 0.5^[p^2] = exp[p^2 log(0.5)].

Spoiler

Now your integral has the form P = FG&integral;_{α}^{β} exp[p^2 log(0.5)] dp. I'm afraid it appears to have no simple antiderivative or closed form solution; the list of integrals on Wikipedia says the closed-form answer would be complex: FGsqrt[pi/(-4olog(0.5))]erf(sqrt[plog(0.5)]), where "erf" is the error function (look it up on Wikipedia).

Yeah, I got as far as step 3 and hit the "Erf" point (this is also, ironically, where I started with the Wolfram Integrator), but I don't really understand how to apply the error function. It seems that it is, itself, an integral that doesn't have a formulaic antiderivative (I'm probably butchering terms).

Edit:Spoiler

Program a routine to do it numerically that takes your bounds and a fineness measure as inputs.

I'm beginning to think that may be the only solution.

[Simon_Jester]

What happens if we define a variable x such that

x = ((r-O)/F)^2

I can tell you right away that won't work. Can you figure out how to quickly see why?

Ah, no. I'm not actually that good at integrals; it was just the only suggestion I could think of in a reasonable time frame. Thank you for explaining why it wouldn't work.

[Steel]

I don't really understand how to apply the error function. It seems that it is, itself, an integral that doesn't have a formulaic antiderivative

If the answer had been Cos or sign you would be happy right?

Well how the hell do you calculate them without either using a calculator or knowing half a dozen values? You could very well define Sin(x) = Integral[Cos(x)] for the appropriate limits. If you can understand the shape of the graph of cos and sin you can use the analytic result, but you still use something to work out specific values, and thats going to be the same here. Erf has a simple graph so you can sketch out the solution the same as if it was another function you were familiar with but cant evaluate in your head for the vast majority of values (Such as Sin, Tan, Exp).

Theres nothing wrong with having the definition for a function as f(x) = Integral[something I cant do].

All that matters is you can write it down and have some way to evaluate it. Also putting it in this form will allow you to work out some identities and possibly manipulate similar results. Erf is a common(ish) function and will be in lots of maths programs as it comes up in a lot of contexts, so having got your solution in this form you can get something to graph it easily and as erf is not a complicated function you can actually visualise (roughly) how changing the bits in it will affect the result.

[Kuroneko]

Exponentials in variable-squared are problamatic; that's what the error function is for:
erf(x) = A Int_0^x exp(-t²) dt, A = 2/sqrt(π).
The error function does not have an elementary form. In general, there's no way to calculate it except numerically.

Anyway, first consider an integral of 2^{-r²}. You can change it into a natural exponential with r' = sqrt[log(2)] r, so we should try s = sqrt[log(2)] (r-O)/F. This gives: Spoiler

Int 2^{-((r-O)/F)²} dr = B erf(s), s = sqrt[log(2)](r-O)/F, B = (F/2)sqrt(π/log 2)

This is purely real, not complex.

[Darth Yoshi]

Well, substituting (r-O)/F with x gives an intregrand of F*0.5x^2. Substituting again with u=x2 gives F(0.5u/2x). Couldn't you evaluate in the form of udx then?

P=GF(0.5u+1ln(x) - integral((u+1)0.5uln(x)du)
P=GF(0.5u+1ln(x) - integral((u+1)0.5uln(√u)du)

Or would that eventually lead to a dead end?

[sketerpot]

Apparently integrals.com can only handle indefinite integrals. According to Wolfram Alpha, making that a definite integral doesn't help you much. Here's the solution it comes up with (which you need to multiply by G, by the way):



And if your calculator can't handle erf, use a computer that can.

[McC]

Out of curiosity, will a computer calculate erf with a substantially higher degree of precision than calculating the integral using very small increments of a and b following the trapezoid rule?

Bearing in mind that r, O, and F tend to go into the tens or hundreds of thousands and I would (well, I have, at this point) use a and b within increments of 1 (i.e. (0,1),(1,2),(2,3), etc.).

[Surlethe]

Any sort of casual integration program will probably use the trapezoid rule, or some variant of it, just with many more intervals than you will. Whether that will substantially increase the precision of the calculation, I couldn't say, but simply programming a computer to do the calculation for you will save a lot of time.

Also, if you're dealing with extreme cases of p = (r - O)/F (i.e., r dominates O and F), you might be able to get by with approximating erf by 1 or -1.

[Kuroneko]

Making a fixed step size h~1 is not a good idea. If you do that, you would almost certainly be better off just using erf from some math package (it's standard in C, for example). If for some reason you have to calculate it yourself, consider recursively adding more points and check if there's any actual improvement going on.

I'll illustrate what I mean using Simpson's rule, although you it's even simpler for trapezoidal rule. If you recall the meat of Simpson's rule,
[1] z = f(a) + 4(a+1h) + 2f(a+2h) + 4f(a+3h) + ... + 4f(a+(n-1)h) + f(b),
and ask yourself: what happens when you triple the number of subintervals, using H = h/3 instead? Then of course you would just use
[2] Z = f(a) + 4(a+1H) + 2f(a+2H) + 4f(a+3H) + ... {4,2}f(a+kH) + ... + 4f(a+(N-1)H)) + f(b),
but notice that a+3jH = a+jh, and they already have the correct coefficients, meaning you don't even have to save your values--all you have to do is skip the k = 0 mod 3 points when calculating Z, and just add z to it. This assumes n was divisible by 3 (and since Simpson's rule requires n to be even, you must start with n a multiple of 6).

So the numerical integral based on the z-value would be w = (hz)/3, while on the Z-value W = (HZ)/3. The recursion would stop when they're close together in either the absolute (|W-w|<ε) or relative (|1-W/w|<ε) sense, or both, i.e., when using more points doesn't seem to improve your accuracy very much, or when some maximum number of recursion steps has been reached.

For trapezoidal rule, you can do the same thing, doubling the number of subintervals, and ignoring those with k = 0 mod 2 on the next iteration (start with n even). It's a pretty good way to automatically detect how good your numerical estimate is without making too many unnecessary function calls.

Most numerical integration packages are even cleverer than this, though, by trying to automatically detect which subintervals are the most problematic. But in your case, you already know which: the most interesting things happen around O.