List Archives >  Maple User Group List Archive >  Archive by date >  This Month By Date >  This Month By Topic

[MUG] Integration problem

Search email archive for  

[MUG] Integration problem
Author: Ben Whitaker    Posted: 26/03/2001 19:30:37 GDT
>> From: "Ben Whitaker" "b.j.whitaker" Hi, I am trying to solve a problem that crops up in simple molecular orbital theory. I am interested in trying to obtain an expression for the overlap integral between two hydrogen 1s atomic orbitals. I have pasted in the worksheet that I am developing at the bottom of this message but in text it is almost impossible to read so I don't know how helpful it will be. My problem is that Maple gets the integration wrong. In polar coordinates one ends up with an integral that has a discontinuity in it so the obvious thing to do is to split it into two parts. One ends up with two integrands exp(-ra - R) (ra + R + 1) + exp(ra - R) (R - ra + 1) A1 := ---------------------------------------------------- R ra and exp(-ra + R) (ra - R + 1) + exp(-ra - R) (ra + R + 1) A2 := ----------------------------------------------------- R ra A1 if ra is less than R; A2 if ra is greater than R which you need to evaulate for ra from 0 to infinity. Integrating A1 from 0 to R and A2 from R to infinity and adding the results using Maple gives the wrong answer. Maple (version 5 and 6.02) gets 3 2 exp(-R) (R + 3 R + 6 R + 6) 1/3 ----------------------------- R but the correct answer is 2 1/3 exp(-R) (R + 3 R + 3) (which you can get by hand or doing the integral in elliptic coordinates with Maple) My question is: Is there a trick that I don't know about that I should use when dealing with the integration of discontinuous functions and why does Maple get the answer wrong using the first method. Regards Ben Whitaker *********************************************************************** Dr Benjamin J Whitaker tel: +44 113 233 6580 School of Chemistry fax: +44 113 233 6565 University of Leeds web: www.chem.leeds.ac.uk Leeds LS2 9JT, UK *********************************************************************** Here for the interested is a text version of the worksheet > restart; We need to evaluate S = int(psi(1*Sa)*psi(1*Sb),r) but as we have already seen the problem is that the two functions have different origins. Let's put the origin on atom A then the integral we need to do is Pi*int(exp(-r[a])*exp(-r[b]),r[a]) but we need to remember that the volume element dra expressed in spherical co-ordinates is dra = r^2*sin(theta)*dra*dtheta*dphi so > S:=(1/Pi)*Int((ra)^2*exp(-ra),ra=0..infinity)*Int(1,phi=0..2*Pi)*Int(sin(the ta)*exp(-rb),theta=0..Pi); We can use the cosine rule to express rb in terms of ra > rb:=sqrt(ra^2+R^2-2*R*ra*cos(theta)); So the overlap integral now becomes > S; infinity 2 Pi Pi / / / | 2 | | | ra exp(-ra) dra | 1 dphi | | | | / / / 0 0 0 2 2 sin(theta) exp(-sqrt(ra + R - 2 R ra cos(theta))) dtheta/Pi Let's now look at the integral with respect to theta. We have > A:=int(sin(theta)*exp(-rb),theta=0..Pi); 2 2 A := - (exp(-sqrt((ra + R) )) sqrt((ra + R) ) 2 + exp(-sqrt((ra + R) )) 2 2 - exp(-sqrt((ra - R) )) sqrt((ra - R) ) 2 - exp(-sqrt((ra - R) )))/(R ra) We have a bit of a problem here because the value of this integral depends on whether ra is bigger or smaller than R. We could probably get Maple to work this out for us but it easier to do the next bit by inspection. If ra< R then > A1:=(1/(R*ra))*(exp(-(R+ra))*(ra+R+1)+exp(-(R-ra))*(R-ra+1)); exp(-ra - R) (ra + R + 1) + exp(ra - R) (R - ra + 1) A1 := ---------------------------------------------------- R ra or if ra > R we have > A2:=(1/(R*ra))*(exp(-(ra-R))*(ra-R+1)+exp(-(R+ra))*(R+1+ra)); exp(-ra + R) (ra - R + 1) + exp(-ra - R) (ra + R + 1) A2 := ----------------------------------------------------- R ra Now S = int(ra^2*exp(-ra),ra = 0 .. infinity)*int(A,phi = 0 .. 2*Pi) where A = A1 for ra = 0 to R and A2 for ra = R to infinity. Since the integrand A doesn't depend on phi the integral over phi just gives us a factor of 2*pi. The factor Pi will cancel with the 1/Pi normalization factor in the original integral so we are left with following two halves for the overal integral > S1:=int(ra^2*exp(-ra)*2*A1,ra=0..R); S2:=int(ra^2*exp(-ra)*2*A2,ra=R..infinity); 2 S1 := 1/6 (-15 exp(-3 R) R - 6 exp(-3 R) - 12 exp(-3 R) R 3 2 + 2 exp(-3 R) R exp(2 R) + 6 exp(-3 R) exp(2 R) R + 3 exp(-R) R + 6 exp(-R))/R 2 (2 exp(2 R) + 5 R + 2 + 3 exp(2 R) R + 4 R ) exp(-3 R) S2 := 1/2 ------------------------------------------------------ R > simplify(S1+S2); 3 2 exp(-R) (R + 3 R + 6 R + 6) 1/3 ----------------------------- R This is not the right answer! What went wrong here? The answer should be S(R) = exp(-R)*(1+R+R^2/3) . Maple has a probem with the integration brcause of the discontinuity at ra = R. If we do the integral by hand we can see this. Start with expressions for the two halves of the intergral. 2*Int(ra^2*exp(-ra)*(exp(-(R+ra))*(ra+R+1)+exp(-(R-ra))*(R-ra+1))/(R*ra),ra = 0 .. R)+2*Int(ra^2*exp(-ra)*(exp(-(ra-R))*(ra-R+1)+exp(-(R+ra))*(R+1+ra))/(R*ra), ra = R .. infinity); = 2*exp(-R)*Int(ra*(R+1-ra),ra = 0 .. R)/R+2*exp(-R)*Int(exp(-2*ra)*ra*(ra-R+1),ra = R .. infinity)/R; = exp(-R)*(R^2/3+R)-2*exp(R)*(-R^2*exp(-2*R)/2+exp(-2*R)*(-2*R-1)/4-R*exp(-2*R )*(-2*R-1)/4+exp(-2*R)*(-2*R-1)/4)/R-exp(-R)*(1/2+1/R); = exp(-R)*(R^2/3+R+R+1+1/(2*R)-R-1/2+1+1/(2*R)-1/2-1/R); = exp(-R)*(1+R+R^2/3) Which is what we wanted but what a nightmare! There is an easier way to do this problem. Two-centre integrals using elliptic coordinates > restart: If we work in elliptic coordinates, a point in space is given by the three coordinates > lambda = (ra+rb)/R; > mu = (ra-rb)/R; and phi; the angle that the (ra,rb, R) triangle makes about the interfocal axis. You might want to sketch the coordinate system. Lines of constant mu are ellipses around the internuclaer axis, whie lines of constant lambda are curved lines that cut through the lines of constant mu at 90 degrees. (Think about using Maple's geometry package to help you). The differential volume element in elliptic coordinates is dr := R^3*(lambda^2-mu^2)*d*lambda*d*mu*d*phi/8 We now evaluate the overlap integral int(exp(-r[a])*exp(-r[b]),r)/Pi Given the definitions of lambda and mu it should be clear that since ra+rb can never be less than R we must have 1 <= lambda and lambda < infinity. Similarly ra-rb; can never be greater than R so -1 <= mu and mu <= 1. From the definition of lambda it follows that ra+rb = R*lambda, so the integral we need to evaluate expressed in elliptic cordinates is > S:=(1/Pi)*Int(Int(Int((R^3/8)*(lambda^2-mu^2)*exp(-R*lambda),mu=-1..1),lambd a=1..infinity),phi=0..2*Pi); The integral is a lot simpler than before but even now Maple gets indigestion if we attempt it all at once, so let's do it in bits > Imu:=int((R^3/8)*(lambda^2-mu^2)*exp(-R*lambda),mu=-1..1); 3 3 2 Imu := - 1/12 R exp(-R lambda) + 1/4 R lambda exp(-R lambda) > factor(Imu); 3 2 1/12 R exp(-R lambda) (-1 + 3 lambda ) > Ilambda:=int(Imu,lambda=1..infinity); Definite integration: Can't determine if the integral is convergent. Need to know the sign of --> R Will now try indefinite integration and then take limits. 2 Ilambda := lim 1/12 R exp(-R lambda) lambda -> infinity 2 2 - 1/4 R lambda exp(-R lambda) - 1/2 R lambda exp(-R lambda) - 1/2 exp(-R lambda) 2 + 1/6 R exp(-R) + 1/2 R exp(-R) + 1/2 exp(-R) So we see where Maple has the problem. It's easy to fix: > assume(R>0): > S:=int((1/Pi)*eval(Ilambda),phi=0..2*Pi); 2 S := 2 (1/6 R~ + 1/2 R~ + 1/2) exp(-R~) > simplify(S); 2 1/3 exp(-R~) (R~ + 3 R~ + 3) Whch is the answer we wanted. *********************************************************************** Dr Benjamin J Whitaker tel: +44 113 233 6580 School of Chemistry fax: +44 113 233 6565 University of Leeds web: www.chem.leeds.ac.uk Leeds LS2 9JT, UK ***********************************************************************

[MUG] Integration problem
Author: Ben Whitaker    Posted: 29/03/2001 19:30:37 GDT
>> From: "Ben Whitaker" "b.j.whitaker" >> From: "Ben Whitaker" "b.j.whitaker" Hi, I am trying to solve a problem that crops up in simple molecular orbital theory. I am interested in trying to obtain an expression for the overlap integral between two hydrogen 1s atomic orbitals. I have pasted in the worksheet that I am developing at the bottom of this message but in text it is almost impossible to read so I don't know how helpful it will be. My problem is that Maple gets the integration wrong. In polar coordinates one ends up with an integral that has a discontinuity in it so the obvious thing to do is to split it into two parts. One ends up with two integrands exp(-ra - R) (ra + R + 1) + exp(ra - R) (R - ra + 1) A1 := ---------------------------------------------------- R ra and exp(-ra + R) (ra - R + 1) + exp(-ra - R) (ra + R + 1) A2 := ----------------------------------------------------- R ra A1 if ra is less than R; A2 if ra is greater than R which you need to evaulate for ra from 0 to infinity. Integrating A1 from 0 to R and A2 from R to infinity and adding the results using Maple gives the wrong answer. Maple (version 5 and 6.02) gets 3 2 exp(-R) (R + 3 R + 6 R + 6) 1/3 ----------------------------- R but the correct answer is 2 1/3 exp(-R) (R + 3 R + 3) (which you can get by hand or doing the integral in elliptic coordinates with Maple) My question is: Is there a trick that I don't know about that I should use when dealing with the integration of discontinuous functions and why does Maple get the answer wrong using the first method. Regards Ben Whitaker *********************************************************************** Dr Benjamin J Whitaker tel: +44 113 233 6580 School of Chemistry fax: +44 113 233 6565 University of Leeds web: www.chem.leeds.ac.uk Leeds LS2 9JT, UK *********************************************************************** Here for the interested is a text version of the worksheet > restart; We need to evaluate S = int(psi(1*Sa)*psi(1*Sb),r) but as we have already seen the problem is that the two functions have different origins. Let's put the origin on atom A then the integral we need to do is Pi*int(exp(-r[a])*exp(-r[b]),r[a]) but we need to remember that the volume element dra expressed in spherical co-ordinates is dra = r^2*sin(theta)*dra*dtheta*dphi so > S:=(1/Pi)*Int((ra)^2*exp(-ra),ra=0..infinity)*Int(1,phi=0..2*Pi)*Int(sin(the ta)*exp(-rb),theta=0..Pi); We can use the cosine rule to express rb in terms of ra > rb:=sqrt(ra^2+R^2-2*R*ra*cos(theta)); So the overlap integral now becomes > S; infinity 2 Pi Pi / / / | 2 | | | ra exp(-ra) dra | 1 dphi | | | | / / / 0 0 0 2 2 sin(theta) exp(-sqrt(ra + R - 2 R ra cos(theta))) dtheta/Pi Let's now look at the integral with respect to theta. We have > A:=int(sin(theta)*exp(-rb),theta=0..Pi); 2 2 A := - (exp(-sqrt((ra + R) )) sqrt((ra + R) ) 2 + exp(-sqrt((ra + R) )) 2 2 - exp(-sqrt((ra - R) )) sqrt((ra - R) ) 2 - exp(-sqrt((ra - R) )))/(R ra) We have a bit of a problem here because the value of this integral depends on whether ra is bigger or smaller than R. We could probably get Maple to work this out for us but it easier to do the next bit by inspection. If ra< R then > A1:=(1/(R*ra))*(exp(-(R+ra))*(ra+R+1)+exp(-(R-ra))*(R-ra+1)); exp(-ra - R) (ra + R + 1) + exp(ra - R) (R - ra + 1) A1 := ---------------------------------------------------- R ra or if ra > R we have > A2:=(1/(R*ra))*(exp(-(ra-R))*(ra-R+1)+exp(-(R+ra))*(R+1+ra)); exp(-ra + R) (ra - R + 1) + exp(-ra - R) (ra + R + 1) A2 := ----------------------------------------------------- R ra Now S = int(ra^2*exp(-ra),ra = 0 .. infinity)*int(A,phi = 0 .. 2*Pi) where A = A1 for ra = 0 to R and A2 for ra = R to infinity. Since the integrand A doesn't depend on phi the integral over phi just gives us a factor of 2*pi. The factor Pi will cancel with the 1/Pi normalization factor in the original integral so we are left with following two halves for the overal integral > S1:=int(ra^2*exp(-ra)*2*A1,ra=0..R); S2:=int(ra^2*exp(-ra)*2*A2,ra=R..infinity); 2 S1 := 1/6 (-15 exp(-3 R) R - 6 exp(-3 R) - 12 exp(-3 R) R 3 2 + 2 exp(-3 R) R exp(2 R) + 6 exp(-3 R) exp(2 R) R + 3 exp(-R) R + 6 exp(-R))/R 2 (2 exp(2 R) + 5 R + 2 + 3 exp(2 R) R + 4 R ) exp(-3 R) S2 := 1/2 ------------------------------------------------------ R > simplify(S1+S2); 3 2 exp(-R) (R + 3 R + 6 R + 6) 1/3 ----------------------------- R This is not the right answer! What went wrong here? The answer should be S(R) = exp(-R)*(1+R+R^2/3) . Maple has a probem with the integration brcause of the discontinuity at ra = R. If we do the integral by hand we can see this. Start with expressions for the two halves of the intergral. 2*Int(ra^2*exp(-ra)*(exp(-(R+ra))*(ra+R+1)+exp(-(R-ra))*(R-ra+1))/(R*ra),ra = 0 .. R)+2*Int(ra^2*exp(-ra)*(exp(-(ra-R))*(ra-R+1)+exp(-(R+ra))*(R+1+ra))/(R*ra), ra = R .. infinity); = 2*exp(-R)*Int(ra*(R+1-ra),ra = 0 .. R)/R+2*exp(-R)*Int(exp(-2*ra)*ra*(ra-R+1),ra = R .. infinity)/R; = exp(-R)*(R^2/3+R)-2*exp(R)*(-R^2*exp(-2*R)/2+exp(-2*R)*(-2*R-1)/4-R*exp(-2*R )*(-2*R-1)/4+exp(-2*R)*(-2*R-1)/4)/R-exp(-R)*(1/2+1/R); = exp(-R)*(R^2/3+R+R+1+1/(2*R)-R-1/2+1+1/(2*R)-1/2-1/R); = exp(-R)*(1+R+R^2/3) Which is what we wanted but what a nightmare! There is an easier way to do this problem. Two-centre integrals using elliptic coordinates > restart: If we work in elliptic coordinates, a point in space is given by the three coordinates > lambda = (ra+rb)/R; > mu = (ra-rb)/R; and phi; the angle that the (ra,rb, R) triangle makes about the interfocal axis. You might want to sketch the coordinate system. Lines of constant mu are ellipses around the internuclaer axis, whie lines of constant lambda are curved lines that cut through the lines of constant mu at 90 degrees. (Think about using Maple's geometry package to help you). The differential volume element in elliptic coordinates is dr := R^3*(lambda^2-mu^2)*d*lambda*d*mu*d*phi/8 We now evaluate the overlap integral int(exp(-r[a])*exp(-r[b]),r)/Pi Given the definitions of lambda and mu it should be clear that since ra+rb can never be less than R we must have 1 <= lambda and lambda < infinity. Similarly ra-rb; can never be greater than R so -1 <= mu and mu <= 1. From the definition of lambda it follows that ra+rb = R*lambda, so the integral we need to evaluate expressed in elliptic cordinates is > S:=(1/Pi)*Int(Int(Int((R^3/8)*(lambda^2-mu^2)*exp(-R*lambda),mu=-1..1),lambd a=1..infinity),phi=0..2*Pi); The integral is a lot simpler than before but even now Maple gets indigestion if we attempt it all at once, so let's do it in bits > Imu:=int((R^3/8)*(lambda^2-mu^2)*exp(-R*lambda),mu=-1..1); 3 3 2 Imu := - 1/12 R exp(-R lambda) + 1/4 R lambda exp(-R lambda) > factor(Imu); 3 2 1/12 R exp(-R lambda) (-1 + 3 lambda ) > Ilambda:=int(Imu,lambda=1..infinity); Definite integration: Can't determine if the integral is convergent. Need to know the sign of --> R Will now try indefinite integration and then take limits. 2 Ilambda := lim 1/12 R exp(-R lambda) lambda -> infinity 2 2 - 1/4 R lambda exp(-R lambda) - 1/2 R lambda exp(-R lambda) - 1/2 exp(-R lambda) 2 + 1/6 R exp(-R) + 1/2 R exp(-R) + 1/2 exp(-R) So we see where Maple has the problem. It's easy to fix: > assume(R>0): > S:=int((1/Pi)*eval(Ilambda),phi=0..2*Pi); 2 S := 2 (1/6 R~ + 1/2 R~ + 1/2) exp(-R~) > simplify(S); 2 1/3 exp(-R~) (R~ + 3 R~ + 3) Whch is the answer we wanted. *********************************************************************** Dr Benjamin J Whitaker tel: +44 113 233 6580 School of Chemistry fax: +44 113 233 6565 University of Leeds web: www.chem.leeds.ac.uk Leeds LS2 9JT, UK ***********************************************************************

[MUG] Re: Integration problem
Author: Maple Group    Posted: 30/03/2001 22:51:42 GDT
>> From: Maple Group "maple_gr" -=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=- Date: Wed, 28 Mar 2001 13:59:36 +0200 From: Helmut Kahovec "helmut.kahovec" To: "maple-list" Subject: Integration problem You have to carefully use the assume facility. This is Maple6 (I skipped some of the intermediate results): > restart; > rb:=sqrt(ra^2+R^2-2*R*ra*cos(theta)); > S:=``(1/Pi)*Int(1,phi=0..2*Pi)* > Int( > Int(sin(theta)*exp(-rb),theta=0..Pi)*(ra)^2*exp(-ra), > ra=0..infinity > ); > R:=R1: # <===!!!=== > assume(ra<R1); # <===!!!=== > Int(sin(theta)*exp(-rb),theta=0..Pi); > A1:=simplify(value(%),assume=positive); > Int(A1*(ra)^2*exp(-ra),ra=0..R1); > B1:=simplify(value(%),assume=positive); > S1:=2*B1; S1 := 1/6 exp(-R1) (15 exp(-2 R1) R1 + 6 exp(-2 R1) 2 3 2 + 12 exp(-2 R1) R1 + 2 R1 + 6 R1 - 3 R1 - 6)/R1 > ra,R:='ra','R': # <===!!!=== > R:=R2: # <===!!!=== > assume(ra>R2); # <===!!!=== > Int(sin(theta)*exp(-rb),theta=0..Pi); > A2:=simplify(value(%),assume=positive); > Int(A2*(ra)^2*exp(-ra),ra=R2..infinity); > B2:=simplify(value(%),assume=positive); > S2:=2*B2; S2 := 2 -((5 R2 + 2 + 4 R2 - 3 R2 exp(2 R2) - 2 exp(2 R2)) exp(-3 R2))/R2/2 > R:='R': # <===!!!=== > simplify(eval(S1+S2,{R1=R,R2=R})); # <===!!!=== 2 1/3 exp(-R) (R + 3 R + 3) With kind regards, Helmut -=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=- From: Samn Sherwood W Civ AFRL/HEDB "Sherwood.Samn" To: "''" "maple-list" Subject: Integration problem Date: Wed, 28 Mar 2001 08:02:43 -0600 Dear Ben, A1, as it is written, diverges. Also, do you mean integrating A2 from R to infinity? Sherwood Samn

[MUG] Integration problem
Author: Benj FitzPatrick    Posted: Tue, 3 Dec 2002 13:16:10 -0800 (
>> From: Benj FitzPatrick "benjfitz" I am trying to numerically integrate a certain function, say f(x,y). It turns out that x and y are functions of t, so the overall integration is done with respect to t. Is there any way to do the integration and place constraints on x or y (I want to be able to have a max value for x that is used when x exceeds the limit). If more particulars would be helpful let me know. __________________________________________________ Do you Yahoo!? Yahoo! Mail Plus - Powerful. Affordable. Sign up now. http://mailplus.yahoo.com

[MUG] Re: Integration problem
Author: Maple User Group    Posted: Fri, 6 Dec 2002 18:00:16 -0500 (
>> From: Maple User Group "maple_gr" On Tue, 3 Dec 2002, Benj FitzPatrick wrote: | I am trying to numerically integrate a certain | function, say f(x,y). It turns out that x and y are | functions of t, so the overall integration is done | with respect to t. Is there any way to do the | integration and place constraints on x or y (I want to | be able to have a max value for x that is used when x | exceeds the limit). -=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=- Date: Wed, 4 Dec 2002 14:58:14 -0500 (EST) From: Carl Devore "devore" To: "maple-list" Subject: Integration problem Use `min`. For example, f:= (x,y)-> x*y: x:= exp(t): y:= GAMMA(t): Int(f(min(1000,x),min(2000,y)), t= 0..100); evalf(%); -=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=- From: "borisal" (Boris Alekseev) To: "maple-list" Date: Thu, 05 Dec 2002 22:24:58 +0300 Subject: Integration problem To insert a constraint on maximal value M of variable x use just min(x,M). f:=(x,y)->x+y; x:=t; y:=t^2; int(f(min(x,M),y),t=a..b); f := (x, y) -> x + y x := t 2 y := t /{ 3 2 \ |{ 1/3 b + 1/2 b b <= M| |{ | |{ 3 2 | \{ 1/3 b + b M - 1/2 M M < b / /{ 3 2 \ |{ - 1/3 a - 1/2 a a <= M| + |{ | |{ 3 2 | \{ - 1/3 a - a M + 1/2 M M < a / In general situation with some region D on the plane (x,y) use characteristic function g(x,y) of region D, which is equal 1 within D and 0 otherwise, to form an integrant f(x,y)*g(x,y) with support on D and invoke int(f(x,y)*g(x,y),t=-infinity..infinity) -=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=-=*=- From: "Dr Francis J. Wright" "F.J.Wright" To: "benjfitz" Subject: Integration problem Date: Fri, 6 Dec 2002 14:18:32 -0000 Can't you just integrate f(min(x,xmax),y) instead of f(x,y)? Francis

Previous by date: [MUG] Syntax question -- updating old script, Kathy Gerber
Next by date: [MUG] Re: Drawing figures (not plots) in Maple, Maple Group
Previous thread: [MUG] Formatting an array for printing, JReinmann524
Next thread: [MUG] Drawing figures (not plots) in Maple, Gaurav Suri