NUMERICAL METHODS/Quantum Mechanics -- Sage

# MATHEMATICAL METHODS OF Science and Engineering # Lectures on Quantum Mechanics #http://www1.uprh.edu/rbaretti #http://www1.uprh.edu/rbaretti/Methodsoftheoreticalphysics.htm # Lectures on Quantum Mechanics http://www1.uprh.edu/rbaretti/LQMIntro.htm # http://www1.uprh.edu/rbaretti/LQMch1.htm #http://www1.uprh.edu/rbaretti/LQMch2.htm # http://www1.uprh.edu/rbaretti/LQMch3.htm # http://www1.uprh.edu/rbaretti/LQMch4.htm # http://www1.uprh.edu/rbaretti/LQMch5.htm # 10 chapters.. #http://www1.uprh.edu/rbaretti/LQMch10.htm

sage: g = Bessel(2); g J_{2} sage: print g J-Bessel function of order 2 sage: g.plot(0,10)

J-Bessel function of order 2

a=2/3 print n(a)

0.666666666666667

print n(pi)

3.14159265358979

var('y,x'); integral(y*exp(-y),y)

-(y + 1)*e^(-y)

var('x,y'); f=(1+x*y)*sin(y); integral(f,y,0,x)

Traceback (click to the left of this block for traceback)
...
Is  x  positive, negative, or zero?

f(x)= x*(sin(x)-x*cos(x))-cos(x) +1 ; diff(f,x)

x |--> x^2*sin(x) - x*cos(x) + 2*sin(x)

var('x'); f2=x-1/2; integral(f2^2,x,0,1)

1/12

phi1=(12)^(1/2)*(x-1/2);integral(phi1^2,x,0,1)

f2=x^2; phi0=1; phi1=(12)^(1/2)*(x-1/2); #integral(f2*phi0,x,0,1) phi2=x^2-(1/3)*phi0-(sqrt(3)/6)*phi1; #phi2=sqrt(540/23)*(x^2-(1/3)*phi0-(1/(6*sqrt(3)))*phi1 ); #integral(phi1*phi2,x,0,1) phi2=sqrt(180)*(x^2-(sqrt(3)/6)*phi1-1/3); integral(phi1*phi2,x,0,1)

a=1/6*sqrt(3); n(a)

0.288675134594813

#data g0,g1,g2 ,t0,t1,t2/8.,14.,8.,0.,50.,100./ var('t');g0=8;g1=14;g2=8;t0=0;t1=50;t2=100; gama=g0*(t-t1)*(t-t2)/((t0-t1)*(t0-t2))+g1*(t-t0)*(t-t2)/((t1-t0)*(t1-t2))+g2*(t-t0)*(t-t1)/((t2-t0)*(t2-t1)); y=plot(gama,t,0,100) show(y)

var('r'); integral(4*pi*exp(-2*r)*r^2,r,0,oo)

pi

s,t=var('s,t') x(s)=exp(-s) integral(x(s)*exp(s*t),t,0,1)

-e^(-s)/s + 1/s

f(x)=(1/x)*(1 - exp(-x) ) ,p0(x) = 1; p1(x)= x; p2(x)=(1/2)*(3*x^2-1)

line 3
SyntaxError: keyword can't be an expression (, line 3)

var('Z'); assume(Z>0); phi(r)=(Z^3/pi)^(1/2)*exp(-Z*r); integral(4*pi*r^4*phi(r)^2,r,0,oo)

3/Z^2

var('r,theta');assume(r>0); #f(theta)=(1/(pi*r))*(1-cos(theta))^(-1/2); f(theta)=1/(pi*r)*(1-cos(theta))^(-1/2); delta=pi/180; integral(f(theta),theta,delta,pi)

integrate(1/sqrt(-cos(theta) + 1), theta, 1/180*pi, pi)/(pi*r)

var('Z,r'); assume (Z>0); v(r)=1/r - ( Z +1/r )*exp(-2*Z*r); phi1s(r)=(Z^3/pi)^(1/2)*exp(-Z*r); phi2s(r)=(Z^3/(32*pi))^(1/2)*(2 - Z*r )*exp(-Z*r/2); integral(phi1s(r)*phi2s(r)*v(r)*4*pi*r^2,r,0,oo)

4096/64827*sqrt(2)*Z

var('Z,r1,r2'); assume (Z>0);assume(r1>0); phi1s(r)=(Z^3/pi)^(1/2)*exp(-Z*r); phi2s(r)=(Z^3/(32*pi))^(1/2)*(2 - Z*r )*exp(-Z*r/2); va=(1/r1)*integral(phi2s(r)^2*4*pi*r^2,r,0,r1) #va vb=integral(phi2s(r)^2*4*pi*r,r,r1,oo) #vb J=integral(phi1s(r1)^2*(va+vb)*4*pi*r1^2,r1,0,oo) J

Traceback (click to the left of this block for traceback)
...
Is  r1  positive, negative, or zero?

-1/8*((Z^4*r1^4 + 4*Z^2*r1^2 + 8*Z*r1 + 8)*e^(-Z*r1)/Z^3 -
8/Z^3)*Z^3/r1

var('Z,r1,r2'); assume (Z>0);assume(r1>0); phi1s(r)=(Z^3/pi)^(1/2)*exp(-Z*r); phi2s(r)=(Z^3/(32*pi))^(1/2)*(2 - Z*r )*exp(-Z*r/2); va=-(1/8)*((Z^4*r^4 + 4*Z^2*r^2 + 8*Z*r + 8)*e^(-Z*r)/Z^3 - 8/Z^3)*Z^3/r ; vb= (1/8)*(Z^3*r^3 - Z^2*r^2 + 2*Z*r + 2)*Z*e^(-Z*r);v=va+vb; J=integral(phi1s(r)*phi2s(r)*v*4*pi*r^2,r,0,oo) J

512/84375*sqrt(2)*Z

n( 512*sqrt(2)/84375)

0.00858165740960029

var('z'); assume(r1>0);assume(Z>0); phi1s=(Z^3/pi)^(1/2)*exp(-Z*r); phi2s=(Z^3/(32*pi))^(1/2)*(2 - Z*r )*exp(-Z*r/2); va=(1/r1)*integral(4*pi*r^2*phi1s,r,0,r1) va #vb(r)= integral(4*pi*r1*phi1s(r1)^2,r1,r,oo); #J1s1s=integral(phi1s(r)^2*(va(r)+vb(r))*4*pi*r^2,r,0,oo) #J1s1s va

Traceback (click to the left of this block for traceback)
...
Is  r1  positive, negative, or zero?

var('z r');assume(z>0); phi2p=(z^5/(96*pi))^(1/2)*r*exp(-z*r/2); integral(phi2p^2*4*pi*r^2,r,0,oo);

var('z r1 r'); assume(z>0); phi1s=(z^3/pi)^(1/2)*exp(-z*r);phi2p=(z^5/(96*pi))^(1/2)*r*exp(-z*r/2); phi2s=(z^3/(32*pi))^(1/2)*(2 - z*r )*exp(-z*r/2); va=(1/r1)*integral(4*pi*r^2*phi2p^2,r,0,r1); #vb= integral(4*pi*r*phi2p^2,r,r1,oo); #phi1s=(z^3/pi)^(1/2)*exp(-z*r1); #integral((va+vb)*4*pi*r1^2,r1,0,oo); va

Traceback (click to the left of this block for traceback)
...
Is  r1  positive, negative, or zero?

var('a'); assume(a>0); integral(exp(-x^2/a),x,0,oo)

1/2*sqrt(pi)*sqrt(a)

n( (1/2)*sqrt(pi))

0.886226925452758

(2/sqrt(pi))*integral(exp(-x^2),x,0,oo)

# exchange term K(1s,2s) var('z ,r1, r'); assume(z>0); assume(r1>0); phi1s=(z^3/pi)^(1/2)*exp(-z*r); phi2p=(z^5/(96*pi))^(1/2)*r*exp(-z*r/2); phi2s=(z^3/(32*pi))^(1/2)*(2 - z*r )*exp(-z*r/2); va=(1/r1)*integral(4*pi*r^2*phi1s*phi2s,r,0,r1) va #vb= integral(4*pi*r*phi1s*phi2s,r,r1,oo); #phi1s=(z^3/pi)^(1/2)*exp(-z*r1);phi2p=(z^5/(96*pi))^(1/2)*r1*exp(-z*r1/2); #phi2s=(z^3/(32*pi))^(1/2)*(2 - z*r1 )*exp(-z*r1/2); #assume(r1>0); #integral(phi1s*phi2s*(va+vb)*4*pi*r1^2,r1,0,oo)

Traceback (click to the left of this block for traceback)
...
Is  r1  positive, negative, or zero?

var('a');Z=2; E = ( a^2 /2 -Z* a) + (1/4) *( a^2 /2 -Z*a) + (17/81)* a -(16/729)* a; y=plot(E,a,1.4,2.2) show(y)

n(112/2187)

0.0512117055326932

var('a');z=3; E =(9/4)* ( a^2 /2 - z* a) +(5/8)*a + 2 *(17/81)*a -(16/729)* a ; y=plot(E,a,2.0,3.1) show(y)

var('a1 a2 r1 r'); #assume(a1>0); #assume(a2>0); assume (r>0); assume(r1>0); phi1s(r)=(a1^3/pi)^(1/2)*exp(-a1*r); phi2s(r)=(a2^3/(32*pi))^(1/2)*(2 - a2*r )*exp(-a2*r/2); va(r1)=(1/r1)*integral(4*pi*r^2*phi2s(r)^2,r,0,r1) vb(r1)= integral(4*pi*r*phi2s(r)^2,r,r1,oo); integral(phi1s(r)^2*(va(r)+vb(r))*4*pi*r^2,r,0,oo)

(8*a1^4*a2 + 20*a1^3*a2^2 + 12*a1^2*a2^3 + 10*a1*a2^4 +
a2^5)*a1^3/(32*a1^7 + 80*a1^6*a2 + 80*a1^5*a2^2 + 40*a1^4*a2^3 +
10*a1^3*a2^4 + a1^2*a2^5)

var(' a2 '); z=3; a1=z-5/16; j1s1s=(5/8)*a1; k1s2s=-(16/729)*a1; j1s2s(a2)=(8*a1**4*a2 + 20*a1**3*a2**2 + 12*a1**2*a2**3 +10*a1*a2**4 +a2**5)*a1**3/(32*a1**7 + 80*a1**6*a2 + 80*a1**5*a2**2 + 40*a1**4*a2**3 +10*a1**3*a2**4 + a1**2*a2**5); e0(a2)=2.*(a1**2/2.-z*a1)+(1./4.)*(a2**2/2.-z*a2); E=e0(a2) +j1s1s+2.*j1s2s(a2)+k1s2s; y=plot(E,a2,.8,2.0); show(y)

var('a'); Z=4; E= (5/2)*( a^2 /2 - Z*a)+ (5/8)*a +4*(17/81)*a +(77/512)*a - 2*(16/729)*a ; y=plot(E,a,2.6,4); show(y)

var('a1 a2 r1 r'); # integral of J(2s,2s) #assume(a2>0); assume (r>0); assume(r1>0); phi1s(r)=(a1^3/pi)^(1/2)*exp(-a1*r); phi2s(r)=(a2^3/(32*pi))^(1/2)*(2 - a2*r )*exp(-a2*r/2); va(r1)=(1/r1)*integral(4*pi*r^2*phi2s(r)^2,r,0,r1) vb(r1)= integral(4*pi*r*phi2s(r)^2,r,r1,oo); integral(phi1s(r)^2*(va(r)+vb(r))*4*pi*r^2,r,0,oo)

Traceback (click to the left of this block for traceback)
...
Is  a2  positive, negative, or zero?

var('a2'); Z=4;a1=Z-5/16; j1s2s(a2)=(8*a1**4*a2 + 20*a1**3*a2**2 + 12*a1**2*a2**3 +10*a1*a2**4 +a2**5)*a1**3/(32*a1**7 + 80*a1**6*a2 + 80*a1**5*a2**2 + 40*a1**4*a2**3 +10*a1**3*a2**4 + a1**2*a2**5); E= 2*( a1^2 /2 - Z*a1)+ 2*(1/4)*( a2^2 /2 - Z*a2)+(5/8)*a1 +4*j1s2s(a2) +(77/512)*a2 - 2*(16/729)*a1 ; y=plot(E,a2,1.0,3.0); show(y)

var('a2'); Z=4;a1=Z-5/16; j1s2s(a2)=(8*a1**4*a2 + 20*a1**3*a2**2 + 12*a1**2*a2**3 +10*a1*a2**4 +a2**5)*a1**3/(32*a1**7 + 80*a1**6*a2 + 80*a1**5*a2**2 + 40*a1**4*a2**3 +10*a1**3*a2**4 + a1**2*a2**5); E= 2*( a1^2 /2 - Z*a1)+ (1)*(1/4)*( a2^2 /2 - Z*a2)+(5/8)*a1 +2*j1s2s(a2) +0.*(77/512)*a2 - 1*(16/729)*a1 ; y=plot(E,a2,2,3); show(y)

var('Z,r'); Z=2; #V=(1/r)*(1-exp(-2*Z*r) ) - Z*exp(-2*Z*r); V=r; y=plot(V,r,1,4); show(y)

Traceback (click to the left of this block for traceback)
...
SyntaxError: Non-ASCII character '\xc2' in file
/tmp/tmpjLDYm8/___code___.py on line 4, but no encoding declared;
see http://www.python.org/peps/pep-0263.html for details

var('r');Z=2; v=(1/r)*(1-exp(-2*Z*r) ) - Z*exp(-2*Z*r); y=plot(v,r,0.0001,5); show(y)

var ('x'); V(x) = -0.0066*x**3 + 0.0967*x**2 - 0.4768*x + 0.975; y=plot(V(x),x,0,5); show(y)

var('a ,z'); assume(z>0); assume(a>0); f1(x)=(4*z^3)^(1/2)*x*exp(-z*x); f2(x)=(8*z/5)^(1/2)*(1-(3*z/4)*x)*exp(-z*x/2); #integral(f1(x)*f2(x),x,0,oo) d2f2(x)=diff(f2(x),x,2); #d2f1(x)=diff(f1(x),x,2); #integral(f1(x)*(-1/2)*d2f1(x),x,0,oo) #integral((-z/x)*f1(x)^2,x,0,oo) integral(f2(x)*(-1/2)*d2f2(x),x,0,oo)

-11/40*z^2

var('a'); assume (a>0); f1(x)=2*x*exp(-x); f2(x)=(1/8^(1/2))*(2-x)*x*exp(-x/2); integral(f2(x)^2,x,0,oo)

k(x)=-(1/2)*diff(f1(x),x,2) -(1/x)*f1(x)

-x*e^(-x)

var ('x'); f1(x) = 2*x*exp(-x); f2(x)=(1/8^(1/2))*(2-x)*x*exp(-x/2) y=plot(f2(x),x,0,16); show(y)

var('x,a'); assume(a>0); #f(x)=exp(-x^2)/(1+x^2); f(x)=exp(-x^4); integral(f(x),x,-1,+1)

-1/2*I*gamma_incomplete(1/4, 1) + 1/2*I*gamma(1/4)

k=9E9 ; sigma=1/k; pia=3.1416 integral(k*sigma*2*pia*r/(1.1^2+r^2)^(1/2),r,0,10)

56.2994706006

n(6.2832*sqrt(101) - 6.2832)

56.8621785026268

n(-(56.29-56.86)/.1)

5.70000000000000

var(' a'); assume(a>0); f(x)=1/(1+a*x^2); integral(f(x),x,-oo,oo);

pi/sqrt(a)

var('a b c'); f(x)=x^3 +2*a*x^2+b*x +c; d2f(x)= diff(f(x),x,2); x=1; d2f(x)

4*a + 6

integrate(exp(-y^2)/(1+y^2)^(1/2),y,-oo,oo)

integrate(e^(-y^2)/sqrt(y^2 + 1), y, -Infinity, +Infinity)

search_doc("definite integral")

Search Documentation: "definite integral"constructions/calculus.html
reference/sage/calculus/calculus.html
reference/sage/calculus/functional.html
reference/sage/functions/piecewise.html
reference/sage/interfaces/maxima.html
reference/sage/misc/functional.html

% fortran !f90 do i=1,5 print*,'float(i)',float(i) enddo

Traceback (click to the left of this block for traceback)
...
SyntaxError: invalid syntax

%fortran C FILE: FIB1.F SUBROUTINE FIB(A,N) C C CALCULATE FIRST N FIBONACCI NUMBERS C INTEGER N REAL*8 A(N) DO I=1,N IF (I.EQ.1) THEN A(I) = 0.0D0 ELSEIF (I.EQ.2) THEN A(I) = 1.0D0 ELSE A(I) = A(I-1) + A(I-2) ENDIF ENDDO END C END FILE FIB1.F

Traceback (click to the left of this block for traceback)
...
OSError: [Errno 2] No such file or directory: 'fortran_module_1.so'

f(x)=x^2*e^(-3*x^3)/(1+x) a=0 b=infinity integrate(f,a,b)

integrate(x^2*e^(-3*x^3)/(x + 1), x, 0, +Infinity)

f(x)=x^2*exp(-x^2)/(1+x^2) a=0 b=5 print integrate(f(x),x,a,b)

integrate(x^2*e^(-x^2)/(x^2 + 1), x, 0, 5)

var('p, a') assume(p>0) assume (a>0) f(r)=exp(-a*r)*sin(p*r) integrate(f(r),r,0,oo)

p/(a^2 + p^2)

#(2.**(3./2.)/pi)/(p**2+1.)**2 f(p)=(2^(3/2)/pi)/(p^2+1)^2 tk=(1/2)*p^2 integrate(4*pi*p^2*tk*f(p)^2,p,0,oo)

1/2

n(5/(4*pi))

0.397887357729738

f(p)=(2^(3/2)/pi)/(p^2+1)^2 tk=(1/2)*p^2 integrate(4*pi*p^2*f(p)^2,p,0,oo)

var('Z,p') assume(Z>0) assume(p>0) f(r)= 4*pi*(Z^3/pi )^(1/2)*(sin(p*r) /p)* exp(-Z*r)*r #integrate(f(r),r,0,oo) fp(p)=8*sqrt(pi)*Z^(5/2)/(Z^4 + 2*Z^2*p^2 + p^4) integrate(4*pi*p^2*fp(p)^2,p,0,oo)

8*pi^3

var('Z,p') assume(Z>0) fp(p)=8*sqrt(pi)*Z^(5/2)/(Z^4 + 2*Z^2*p^2 + p^4) integrate(4*pi*p^2*fp(p)^2,p,0,oo)

8*pi^3

n(8*pi^3)

248.050213442399

var('x,y'); v(x,y)=x^2 + 3*x*y^2; ex(x,y)=-diff(v(x,y),x); x=1; y=2; ex(x,y)

-14

f(p)=exp(-p^2); integral(f(p)^2*4*pi*p^2,p,0,oo)

1/4*pi^(3/2)*sqrt(2)

f(p)= (2/pi)^(3/4)*exp(-p^2) integral(f(p)^2*4*pi*p^2,p,0,oo)

g2(p)=exp(-2*p^2) g1n(p)=(2/pi)^(3/4)*exp(-p^2) integral(g2(p)*g1n(p)*4*pi*p^2,p,0,oo)

1/9*pi^(3/4)*2^(3/4)*sqrt(3)

g2(p)=exp(-2*p^2) g1n(p)=(2/pi)^(3/4)*exp(-p^2) u2(p)=g2(p)-(1/9)*pi^(3/4)*2^(3/4)*sqrt(3)*g1n(p) #integral(u2(p)*g1n(p)*4*pi*p^2,p,0,oo) integral(u2(p)^2*4*pi*p^2,p,0,oo)

1/432*(27*sqrt(2) - 32)*pi^(3/2)*sqrt(2)

n(1/432*(27*sqrt(2) - 32)*pi^(3/2)*sqrt(2))

0.112722112725355

n((1/9)*pi^(3/4)*2^(3/4)*sqrt(3)/sqrt(0.112722112725355) )

2.27482744465886

n(1/sqrt(0.112722112725355 ))

2.97848515558980

u2n(p)=2.97848515558980*exp(-2*p^2) -2.27482744465886*exp(-p^2) integral(u2n(p)^2*4*pi*p^2,p,0,oo)

1/36*((6.65353036655*sqrt(2) + 15.5245197089)*sqrt(3) -
27.1021591018*sqrt(2))*pi^(3/2)*sqrt(2)*sqrt(3)

n((1/36)*((6.65353036655*sqrt(2) + 15.5245197089)*sqrt(3) - 27.1021591018*sqrt(2))*pi^(3/2)*sqrt(2)*sqrt(3))

n(1.621282/2.978485)

0.544331094499385

g1(p)=(2/pi)^3/4 *exp(-p^2 ) g2(p)=2.978485*exp(-2*p^2) -1.621282*exp(-p^2 ) dv(p)=4*pi*p^2 integral(g2(p)^2*dv(p),p,0,10)

var('Z'); phi(p)= (1/(8*pi^3)^(1/2))*8*sqrt(pi)*Z^(5/2)/(Z^4 + 2*Z^2*p^2 + p^4);phi(p)

2*sqrt(2)*Z^(5/2)/((Z^4 + 2*Z^2*p^2 + p^4)*pi)

Z=1; phi(p)= (2^(3/2)/pi)*Z^(5/2)/(p^2+Z^2)^2 integral(phi(p)^2*4*pi*p^2,p,0,oo)

x = var('x') g1 = plot(cos(20*x)*exp(-2*x), 0, 1) g2 = plot(2*exp(-30*x) - exp(-3*x), 0, 1) show(graphics_array([g1, g2],1,1), xmin=0)

x=var('x') g1=(2/pi)^(3/4)*exp(-x^2 ) g2=2.978485*exp(-2*x^2 ) -1.621282*exp(-x^2) phi1= plot((2^(3/2) /pi) /(x^2 + 1 )^2,0,3) phi2= plot( g1 + 8.26E-2 * g2,0,3) show(graphics_array([phi1, phi2],1,2), xmin=0)

Traceback (click to the left of this block for traceback)
...
SyntaxError: invalid syntax

var('p') g1n(p)=(2/pi)^(3/4)*exp(-p^2 ) u2(p)= p*exp(-p^2) dv(p)=4*pi*p^2 g2(p)= 2.11661285743602 *(u2(p)-1.11951513492025*g1n(p)) norm=integral(g2(p)^2*dv(p),p,0,oo);n(norm)

0.999999999999996

var('p') g1n(p)=(2/pi)^(3/4)*exp(-p^2 ) g2n(p)= 2.116613*(p*exp(-p^2 )-1.119515*(2/pi)^(3/4)*exp(-p^2 ) ) dv(p)=4*pi*p^2 norm=integral(g2n(p)*g1n(p)*dv(p),p,0,oo); n(norm)

2.85573950264961e-7

A=matrix([[ -0.423944443 , 0.0170494914 ],[ 0.0170494914 , 0.0841127038 ]]) A.eigenvalues()

[0.0846842114303050, -0.424515950630305]

NUMERICAL METHODS/Quantum Mechanics

1 day ago by reibaretti

Search Documentation: "definite integral"