Integral problems with polar coordinates and variable substitution

[Edison Bias]

H!

I wonder how to solve:

I=\int_{-\infty}^{\infty}e^{-u^2}\frac{1}{1+Cu} du

I have solved:

\int_{-\infty}^{\infty}e^{-u^2}du

which equals

\sqrt{\pi}

and I solved it with polar coordinates and variable substitution.

Thankful for help!

Edison

 

[haruspex]

If C is real, won't the integral diverge in the vicinity of -1/C?

 

[Edison Bias]

C is supposed to be real and yes there is a singularity for u=-1/C but the question is if that ruins the possibillity to calculate the integral. I don't know but if you can calculate with infiny as limits you may be able to calculate with infinity as singularity, or?

I have made two attempts into calculate the above integral with C=1 and these are my results (using integration by parts):

\int u(x)v'(x)dx=u(x)v(x)-\int u'(x)v(x)dx

so if

v'(x)=\frac{1}{1+u}

and

u(x)=e^{-u^2}

then

\int e^{-u^2} \frac{1}{1+u}du =[e^{-u^2} ln(1+u)]+\int 2u e^{-u^2} ln (1+u)du

or if

v'(x)=e^{-u^2}

and

u(x)=\frac{1}{1+u}

then

\int e^{-u^2} \frac{1}{1+u}du =[F(e^{-u^2})\frac{1}{1+u}]+\int F(e^{-u^2})\frac{1}{(1+u)^2}du

which both are rubbish, right? At least to the extent that integration has not become a single bit more simple. F is by the way the primitive function to the expression within parenthesis. We also have the problem that the first formula brackets will contain a negative ln number which does not excist. However, maybe e^(-u^2) is "faster" :) Otherwise it looks like both bracket-expressions can be omitted due to zero, but I don't know.

Edison

 

[Charles Link]

I think the complex variable residue theorem will get you to the result. Standard calculus substitutions might not work for this integral. As haruspex mentioned though, it might diverge. I see the "pole" in this case is right on the x-axis. That complicates the calculation using the residue theorem.

 

[haruspex]

the question is if that ruins the possibillity to calculate the integral. I don't know but if you can calculate with infiny as limits you may be able to calculate with infinity as singularity, or?

The two are independent. Whether you can integrate through a singularity depends on how fast it runs away to infinity there.
Consider a small region around the singularity. The numerator won't vary much through that, so you can bound it with a couple of constants.

 

[Edison Bias]

Things have changed, now I only need to solve this integral where the singularity problem is gone:

I=\int_{0}^{\infty}e^{-u^2}\frac{1}{1+Cu} du

Edison

 

[Edison Bias]

I have now tried to solve the integral (C=1) numerically using this approximation

I\approx \sum_{k=0}^{N}f(kd)d=\sum_{k=0}^{N}e^{-(kd)^2}\frac{1}{1+kd}d

I have then put d=0,2 och N=10 and this while the last term (e.g at kd=2 to emphasize) becomes as small as e^-(2^2)/(1+2)*0,2=0,13%.

This then gives

f(k=0)d=0,20
f(k=1)d=0,16
f(k=2)d=0,12
f(k=3)d=0,088
f(k=4)d=0,059
f(k=5)d=0,037
f(k=6)d=0,022
f(k=7)d=0,012
f(k=8)d=0,0060
f(k=9)d=0,0029
f(k=10)d=0,0013

and summed up it reaches 0,7082 which is close to 1/sqrt(2).

So it seems that you don't have to get so far from zero until there is nothing more to add. Infinity as limits is therefore academic, or?

Edison

 

[haruspex]

I have now tried to solve the integral (C=1) numerically using this approximation

I\approx \sum_{k=0}^{N}f(kd)d=\sum_{k=0}^{N}e^{-(kd)^2}\frac{1}{1+kd}d

I have then put d=0,2 och N=10 and this while the last term (e.g at kd=2 to emphasize) becomes as small as e^-(2^2)/(1+2)*0,2=0,13%.

This then gives

f(k=0)d=0,20
f(k=1)d=0,16
f(k=2)d=0,12
f(k=3)d=0,088
f(k=4)d=0,059
f(k=5)d=0,037
f(k=6)d=0,022
f(k=7)d=0,012
f(k=8)d=0,0060
f(k=9)d=0,0029
f(k=10)d=0,0013

and summed up it reaches 0,7082 which is close to 1/sqrt(2).

So it seems that you don't have to get so far from zero until there is nothing more to add. Infinity as limits is therefore academic, or?

Edison

Yes, it should converge very quickly. But I doubt the proximity to 1/√2 is significant, given that the d value you used was 0.2. If a d value of 0.5, say, had given that result it would have been more suggestive.

 

[Edison Bias]

I don't know if I agree with you because you can't have too large steps (d) because then the partial sum is not representative, but I do agree with you that using my numerical integration beyond kd=2 is better and your suggestion could be done by putting N=25 instead, at the expense of making the calculations by hand 25 times... (I have a HP48SX from school days but it has been laying around collecting dust for 20 years and I have recently swiched the leaking batteries and begun using it. Someday I will sit down and go through both of the thick instruction manuals to learn how to really use it, again).

Edison

 

[haruspex]

I don't know if I agree with you because you can't have too large steps (d) because then the partial sum is not representative, but I do agree with you that using my numerical integration beyond kd=2 is better and your suggestion could be done by putting N=25 instead, at the expense of making the calculations by hand 25 times... (I have a HP48SX from school days but it has been laying around collecting dust for 20 years and I have recently swiched the leaking batteries and begun using it. Someday I will sit down and go through both of the thick instruction manuals to learn how to really use it, again).

Edison

I would guess that the greater inaccuracy is in summing the values you have, as though the graph were a step function over that range. The trapezium method would be more accurate.

 

[Edison Bias]

I think it is obvious that it converts (is that the right word?) to around

\frac{1}{\sqrt{2}}

Edison

 

[haruspex]

I think it is obvious that it converts (is that the right word?) to around

\frac{1}{\sqrt{2}}

Edison

I agree the value will be close to 1/√2, but it seems likely that is just a coincidence. If you were to get around that number with, say, C=1 and d=1/4 then it would be more suggestive.

 

[Charles Link]

I've been working this one off and on for C=1, and the best thing I have so far is writing $1/(1+u)$ as a geometric series, but these terms (moments of the Gaussian) are alternating and do not seem to converge very quickly...(editing).. And upon taking a closer look- the geometric series is mathematically unsound because $u$ takes on values larger than 1. I'm still at the drawing board on this one. My previous comments about the complex variable residue theorem also don't appear to work for this case for a couple of reasons. Numerical evaluation (like you have been doing) with very small increments might help to hint at what the correct answer may be.

 

[Edison Bias]

I'm honored that you have put so much work into helping me!

I will try to get to understand my HP48SX and then enter my numerical integration formula above and use:

[C=1]
d=0,1
N=50

I think it will be most simple to calculate one term at the time, somehow store them on the stack and finally just add them all up. If I can't find out how to store them on the stack, I will simply write them down with three significant digits. Yet, I will have to write down 50 values :)

Edison

 

[Charles Link]

My latest attempt to solve it looks too complex, but it was worth a try. Just as for the case with C=0, I tried using the $I^2=I_x*I_y$ approach and converting the double integral into polar coordinate form, this time with the integral covering just the first quadrant. The $(1+r \cos(\theta))$ and $(1+r \sin(\theta))$ terms in the denominator look to be problematic though.

 

[Edison Bias]

But for C=0 you actually can solve it as you have done, I think. Let's say

I^2=\int_{-\infty}^{\infty}e^{-x^2}dx*\int_{-\infty}^{\infty}e^{-y^2}dy=\int \int_{-\infty}^{\infty} e^{-(x^2+y^2)}dxdy

where

x^2+y^2=r^2

then we can solve

I^2=\int_{0}^{2\pi} \int_{0}^{\infty} e^{-r^2}rdr d\phi=2\pi \int_{0}^{\infty}e^{-r^2}rdr=\pi\int_{0}^{\infty} e^{-r^2}d(r^2)=[-\pi e^{-r^2}]_{0}^{\infty}=\pi

I'm not sure here but if we stay in the first quadrant the integral should be

I=\sqrt{\frac{\pi}{2}}

due to only moving half a turn (pi and not 2pi) so to speak.

It is interersting to note that integration from zero to infinity actually gets the volume of half the Maxwellian Distribution, or so I understand it.

But is this really correct for us?

Can we really do this because aren't we interested in the integral in one plane only?

Edison

 

[Charles Link]

Yes, it would be correct for us to do, even for C=1 and have the polar coordinates cover the first quadrant. I thought we might get lucky here, but the polar coordinate integral for C=1 in the first quadrant looks unmanageable. For your numerical integration, the EXCEL spreadsheet is a very simple way to do it. You can easily do 1000 pts. with just a couple lines of code, and sum a column of numbers. ( basically a "Do" loop or "For/Next" Loop). (editing...)... I found an on-line version of EXCEL and did 300 increments from x=0 to x=3. $\int (\exp(-x^2)/(1+x)) \, dx$ gives approximately .61. I will try to double-check the computer spreadsheet, but it does appear the answer is significantly different from 1/sqrt(2). And yes, the best I can tell, the spreadsheet gave me an accurate answer. Additional comment on your post #16: If you stay in the first quadrant, the result is sqrt($\pi$)/2.

 

[Edison Bias]

Extremely interesting!

So you mean that the integral actually is less than my preliminary answer? I thought that going beyond kd=2 while using smaller d and more terms (N) would get a higher answer, not a lower one. This while the integrand asymptotically approaches zero and this for high kd only. But you have proven that the answer actually is lower than mine and this while using kd=3 and very much smaller d, interesting. So the integral with C=1 is closer to 1/sqrt(3) than 1/sqrt(2) (these answers just look nicer :) ).

Edison

 

[Charles Link]

(I just edited post #17 as your latest post came in, so please read my additional comments on your post #16.). This integral is simple enough in form, that it could very well have a simple answer, but the other thing that normally, (not always, but normally), comes with a simple answer is a way to get to the simple answer. So far, the attempts to get a closed form solution have lead to complexities. There could be a simple answer to this one, but it is possible this one doesn't have a simple closed form. I do plan on giving it a little more effort, but it is quite possible that any closed form solution for this one is rather complex. Perhaps other readers might have a solution for it, but I have not succeeded at finding it. And to comment on # 18, I essentially used a "kd" of .01 and took 300 terms. (out to x=3). For larger values of "kd", the summation could contain significant errors so that the answer to your sum could be larger than the actual integral. I do think my summation (.61) is accurate to better than +/- .005. (Incidentally, I used the same spreadsheet to compute $\int exp(-x^2) \, dx$ from zero to infinity (with a simple formula change) and got .89)

 

[haruspex]

I'm honored that you have put so much work into helping me!

I will try to get to understand my HP48SX and then enter my numerical integration formula above and use:

[C=1]
d=0,1
N=50

I think it will be most simple to calculate one term at the time, somehow store them on the stack and finally just add them all up. If I can't find out how to store them on the stack, I will simply write them down with three significant digits. Yet, I will have to write down 50 values :)

Edison

As I wrote, the trapezium method is more accurate. But using a simple summation you can get bounds on the integral. If your values are y₀ to y_n at intervals of dx, then because the function is monotonically decreasing a lower bound for the area from 0 to n dx is dxΣ₁ⁿy_i and an upper bound is dxΣ₀^n-1y_i. You can see this quite easily by drawing bounding boxes on the graph.

 

[Edison Bias]

I'm sorry but that is too academic to me, I would rather use my approach with less d and more N.

Edison

 

[Edison Bias]

(I just edited post #17 as your latest post came in, so please read my additional comments on your post #16.). This integral is simple enough in form, that it could very well have a simple answer, but the other thing that normally, (not always, but normally), comes with a simple answer is a way to get to the simple answer. So far, the attempts to get a closed form solution have lead to complexities. There could be a simple answer to this one, but it is possible this one doesn't have a simple closed form. I do plan on giving it a little more effort, but it is quite possible that any closed form solution for this one is rather complex. Perhaps other readers might have a solution for it, but I have not succeeded at finding it. And to comment on # 18, I essentially used a "kd" of .01 and took 300 terms. (out to x=3). For larger values of "kd", the summation could contain significant errors so that the answer to your sum could be larger than the actual integral. I do think my summation (.61) is accurate to better than +/- .005. (Incidentally, I used the same spreadsheet to compute $\int exp(-x^2) \, dx$ from zero to infinity (with a simple formula change) and got .89)

I don't agree with you that the integral with C=0 should amount to sqrt(pi)/2. I am somewhat confident that it should amount to sqrt(pi/2). Because, look at the I^2-part and then use the polar coordinates which gives you sqrt(pi) when integrating from -infinity to +infinity. But here we are just integrating from zero to infinity. So the angle for the polar coordinates does not revolve a full turn, it stops with pi. So if our I^2 then results in pi/2 then I must result in sqrt(pi/2), right? On the other hand, who has said that a "half turn" is representative for our integral? It would be more intuitive if one could say that if this polar stuff is a kind of volume integration then using half the space would mean sqrt(pi)/2.

Edison

 

[Charles Link]

I don't agree with you that the integral with C=0 should amount to sqrt(pi)/2. I am somewhat confident that it should amount to sqrt(pi/2). Because, look at the I^2-part and then use the polar coordinates which gives you sqrt(pi) when integrating from -infinity to +infinity. But here we are just integrating from zero to infinity. So the angle for the polar coordinates does not revolve a full turn, it stops with pi. So if our I^2 then results in pi/2 then I must result in sqrt(pi/2), right? On the other hand, who has said that a "half turn" is representative for our integral? It would be more intuitive if one could say that if this polar stuff is a kind of volume integration then using half the space would mean sqrt(pi)/2.

Edison

In the I^2 integral, the x has limits of zero and infinity and so does the y. (When this calculation of I^2 starts out, the I^2 is written as the product of two separate integrals, but then the integrals are combined as one double integral.) When these limits are combined, the result is just the first quadrant. All quadrants will get the same answer because of the radial symmetry. Thereby I^2=$\pi/4$ for the limits of zero to infinity (for both x and y). (It takes a little practice to understand the limits in double integrals, but hopefully this one makes sense.)

 

[Edison Bias]

I think I'm closer to understanding now because strictly positive integration limits makes us reside in one quadrant only, right? So if all quadrants (full "volume", my note) means I^2=Pi then one fouth volume would need to equal Pi/4. Thank you for that explanation even though, with your help, really looking at the limits actually made it easy. If my "volume" approach is the least appliceble, that is.

Edison
PS
One thought here, if we use a cartesian double integral and look at the limits, I agree with you that integration will take place in one quadrant only. But what happens when we move to polar coordinates? I think I^2 still should be calculated with r going from zero to infinity but what about phi? If full volume phi is 2Pi, then our must be Pi/2 and this actually gives your result as a factor of 4 lower value of I^2.

 

[Charles Link]

I think I'm closer to understanding now because strictly positive integration limits makes us reside in one quadrant only, right? So if all quadrants (full "volume", my note) means I^2=Pi then one fouth volume would need to equal Pi/4. Thank you for that explanation even though, with your help, really looking at the limits actually made it easy. If my "volume" approach is the least appliceble, that is.

Edison
PS
One thought here, if we use a cartesian double integral and look at the limits, I agree with you that integration will take place in one quadrant only. But what happens when we move to polar coordinates? I think I^2 still should be calculated with r going from zero to infinity but what about phi? If full volume phi is 2Pi, then our must be Pi/2 and this actually gives your result as a factor of 4 lower value of I^2.

Yes, in polar coordinates, the phi limits will be zero and $\pi/2$. Precisely correct.

 

[Edison Bias]

Makes me almost proud

 

[Edison Bias]

You might want to know where these integral problems have come from?

n=\int_{-\infty}^{\infty}f(v)dv=A\int_{-\infty}^{0}e^{-\frac{m(v+v_0)^2/2}{kT}}dv+A\int_{0}^{\infty}e^{-\frac{m(v-v_0)^2/2}{kT}}dv=2A\int_{0}^{\infty}e^{-\frac{m(v-v_0)^2/2}{kT}}dv

which describes the distribution of velocities in matter and therefore the kinetic energy in matter where the kinetic energy might explain why for instance water evaporates at (average) temperatures far blow the boiling point.

I have thus revised the Maxwellian distribution function which states

n=\int_{-\infty}^{\infty}Ae^{\frac{mv^2/2}{kT}}dv

where n is the particle density.

You're not supposed to come up with own ideas here at PF but I can't help critisizing that that distribution function gives an average distribution that is most likely for v=0. Sure, there can be positive and negative velocities (did not understand that at first) but to say that particles in average stands still oversimplifies stuff, I think.

So what I have done is change v to v'-v_0 where v_0 equals average speed, this means that there is a probability for v=0 but that is not the most probable speed and the level of probability depends on weather

mv_0^2=kT

is true or not.

Substituting in the integrals above

u^2=\frac{m(v-v_0)^2/2}{kT}

gives

n=2A\sqrt{\frac{2kT}{m}}\int_{0}^{\infty}e^{-u^2}du-2A\sqrt{\frac{2kT}{m}}\int_{0}^{\infty}e^{-u^2}\frac{v_0}{v_0+u\sqrt{\frac{kT}{m}}}du

where the last integral is our integral with C=1 due to the above condition.

Hope I did not bore you.

Edison

 

[Charles Link]

For the Maxwellian, I can give you some additional information which you might find of interest. You wrote a formula for the $v_x$ distribution where $v_x$ is velocity in the x-direction. This will average to zero. (Because the particle can go in both positive and negative x-directions). Likewise for $v_y$ and $v_z$. These can be converted to spherical coordinates where $v^2=v_{x} ^2 +v_{y}^2+v_{z}^2$ and with spherical symmetry $4 \pi v^2 dv=(dv_x) (dv_y)(dv_z)$. $v$ now represents speed and the v is integrated from 0 to + infinity. Computations can be performed to find the most probable speed, the average speed, and the r.m.s. speed. These turn out to be the well-known results $v_{mp}=(2kT/m)^{1/2}$, $v_{avg}=((8/ \pi)kT/m)^{1/2}$, and $v_{rms}=(3kT/m)^{1/2}$. I will need to study your methodology a little more to see if it is at all useful. It did give a very interesting integral which we so far only have a numerical evaluation of it.

 

[Edison Bias]

Just want to add that my way of seeing it is that particles move and this gives them kinetic energy and thus heat so saying that particles in average stand still has to be wrong. Average speed has to be something we can feel.

 

[haruspex]

Just want to add that my way of seeing it is that particles move and this gives them kinetic energy and thus heat so saying that particles in average stand still has to be wrong. Average speed has to be something we can feel.

Average velocity is not the same as average speed.