MomentumRepresentation -- Sage

# MATHEMATICAL METHODS OF Science and Engineering # Lectures on Quantum Mechanics /Density Functional theory #http://www1.uprh.edu/rbaretti #http://www1.uprh.edu/rbaretti/Methodsoftheoreticalphysics.htm #http://www1.uprh.edu/rbaretti/DensityFunctionalTheory4marzo2010.htm #http://www1.uprh.edu/rbaretti/MomentumspaceIntegrationPart213feb2010.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

var('p') g1n(p)=(2/pi)^(3/4)* exp(-p^2 ) u2(p)=exp(-4*p^2) dv(p)=4*pi*p^2 #g2n(p)=(u2(p)-0.354960391545000*g1n(p))/sqrt(0.120090775836133) g2n(p)=-1.02429503946317*(2/pi)^(3/4)*exp(-p^2)+2.88566010141252*e^(-4*p^2) norm=integral(g2n(p)^2*dv(p),p,0,oo);n(norm)

1.00000000000001

A=matrix([[ -0.483184278 , 0.0192039013 ],[ 0.0192039013 ,0.502940059 ]]) A.eigenvalues()

[0.503313896311354, -0.483558115311354]

var('x') #phi(x)= (0.017*x^3 + 0.304*x^2 - 1.0271*x + 0.9953)/sqrt(44.0169094980943) phi(x)=(1/sqrt(1.22604))*1/(x^2+1)^2 dv(x)=4*pi*x^2 #norm=integral(phi(x)^2*dv(x),x,0,3);n(norm) tk=integral((x^2/2)*phi(x)^2*dv(x),x,0,3);n(tk)

0.453877992191000

f(x)=x/(1+exp(-x^2)) d2x(x)=diff(f(x),x,2) x=1/2 n(d2x(x))

0.723098492199436

A=matrix([[ -0.493270665 ,-0.00618162751 ],[ -0.00618162751 ,-0.0781037211 ]]) A.eigenvalues()

[-0.0780117001638897, -0.493362685936110]

#c1,c2= 1. 0.0147767756 znuc=1 #Fig 2 Plot of Φ1s(p)= c1 *g1 (p) + c2 * g2(p) g1(p)=0.9003252/(p^2+znuc^2)^2 g2(p)=0.291716*(1./(p^2+(znuc/2.)^2)^2 -5.26566*g1(p)) dv(p)=4*pi*p^2 norm=integrate(g2(p)^2*dv(p),p,0,10),n(norm)

znuc=1 c1=1 c2= 0.0147767756 g1(p)=0.9003252/(p^2+znuc^2)^2 g2(p)=0.291716*(1./(p^2+(znuc/2.)^2)^2 -5.26566*g1(p)) phi(p)=c1*g1(p)+c2*g2(p) y=plot(phi(p),p,0,3) show(y)

znuc=1 c1=1 c2= -67.1766129 dv(p)=4*pi*p^2 g1(p)=0.9003252/(p^2+znuc^2)^2 g2(p)=0.291716*(1./(p^2+(znuc/2.)^2)^2 -5.26566*g1(p)) phi2s(p)=1.48845E-2 *(c1*g1(p) + c2*g2(p)) #phi(p)=c1*g1(p)+c2*g2(p) y=plot(g1(p)^2,p,0,1.5) show(y)

#limit of int p=1.5 var('p') znuc=1 phi1s(p)= 0.9003252/(p^2+znuc^2)^2 phi2s(p)= -1.0352*p^4 + 9.6359*p^3 - 23.114*p^2 + 20.444*p - 5.7585 dv(p)=4*pi*p^2 norm=integral_numerical(phi1s(p)*phi2s(p)*dv(p),0,1.5,max-points=500)[0]

Traceback (click to the left of this block for traceback)
...
SyntaxError: keyword can't be an expression

# http://sage.math.washington.edu/home/wdj/teaching/calc2-sage/calc2-sage.pdf #sage: integral(f(x),x,0,1) #1 - cos(1) #sage: RR(integral(f(x),x,0,1)) #0.459697694131860 var('x') f(x)=x^2 area=integral(f(x),x,0,2)

0.864664716763387

A=matrix([[ -0.488956273 , 0.0603011362 ],[ 0.0603011362 , -0.0730859786 ]]) A.eigenvalues()

__main__:1: UserWarning: Using generic algorithm for an inexact
ring, which will probably give incorrect results due to numerical
precision issues.
[-0.0645188106435643, -0.497523440956436]

var('z ,r ,p') assume(z>0) phi2s(r)=sqrt(z^3/(32*pi))*(2-z*r)*exp(-z*r/2) integrate(4*pi*(sin(p*r)/p)*phi2s(r)*r,r,0,oo)

32*(4*p^3*z - p*z^3)*sqrt(pi)*sqrt(2)*z^(3/2)/((64*p^6 + 48*p^4*z^2
+ 12*p^2*z^4 + z^6)*p)

z=1 phi2sp=0.0634936393834597*(32*(4*p^3*z - p*z^3)*sqrt(pi)*sqrt(2)*z^(3/2)/((64*p^6 + 48*p^4*z^2 + 12*p^2*z^4 + z^6)*p)) y=plot(phi2sp,p,1.e-3,3) show(y)

var('z ,r ,p') assume(z>0) dv(p)=4*pi*p^2 phi1s(r)=sqrt(z^3/pi)*exp(-z*r) phi2s(r)=sqrt(z^3/(32*pi))*(2-z*r)*exp(-z*r/2) #integrate(4*pi*(sin(p*r)/p)*phi2s(r)*r,r,0,oo) integrate(phis(r)*r,r,0,oo)

32*(4*p^3*z - p*z^3)*sqrt(pi)*sqrt(2)*z^(3/2)/((64*p^6 + 48*p^4*z^2
+ 12*p^2*z^4 + z^6)*p)

var('z ,p') assume(z>0) dvol(p)=4*pi*p^2 norm2=1/4*sqrt(2)/pi^(3/2) phi1s(p)= (1/pi)*2^(3/2)*z^(5/2)/(p^2 + z^2)^2 phi2s(p)=norm2*32*(4*p^3*z - p*z^3)*(2*pi)^(1/2)*z^(3/2)/((64*p^6 + 48*p^4*z^2 + 12*p^2*z^4 + z^6)*p); norm=integrate(phi2s(p)^2*dvol(p),p,0,oo) norm

n=1 c=137.035999 alfa=1/c e(z)=c^2*(1-(1/2)*(z*alfa/n)^2-((z*alfa)^4/(2*n^4))*(n/(1/2)-3/4 )) #e(z)=-z^2/(2*n^2) plot(e(z)-c^2,z,10,40)

assume(a>0) var('a,x') integral(1/(1+a*x^2),x)

arctan(sqrt(a)*x)/sqrt(a)

arctan(-oo)

-1/2*pi

A=matrix([[ 0.607785285 , -1882.86548 ],[ -1882.86548 , 90284216. ]]) A.eigenvalues()

[9.02842160392669e7, 0.568518377317205]

c=137.035999 16712.5820070147-c^2 #-2066.28301491330 #0.000102880057962364-c^2 9.02842160392669e7-c^2 0.568518377317205-c^2

-18778.2965035507

A=matrix([[ -0.85348E+03 ,-0.15060E+02 ],[ -0.15060E+02 , -0.23129E+03 ]]) A.eigenvectors()

Traceback (click to the left of this block for traceback)
...
AttributeError:
'sage.matrix.matrix_generic_dense.Matrix_generic_dense' object has
no attribute 'eigenvectors'

http://sagemath.org/doc/constructions/linear_algebra.html sage: A = matrix(QQ, [[1,1,0],[0,2,0],[0,0,3]]) sage: A [1 1 0] [0 2 0] [0 0 3] sage: A.eigenvalues() [3, 2, 1] sage: A.eigenvectors_right() [(3, [ (0, 0, 1) ], 1), (2, [ (1, 1, 0) ], 1), (1, [ (1, 0, 0) ], 1)]

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

sage: A = matrix(QQ, [[1,3],[3,4]]) sage: A.eigenvalues() sage: A.eigenvectors_right()

[(-0.8541019662496845?, [(1, -0.618033988749895?)], 1),
(5.854101966249684?, [(1, 1.618033988749895?)], 1)]

A = matrix(QQ, [[ -0.87387E+03, 0.36145E+02 ],[ 0.36145E+02 , -0.20003 ]]) A.eigenvalues() A.eigenvectors_right()

[(-875.3628205804243?, [(1, -0.04130088754805030?)], 1),
(1.292790580424279?, [(1, 24.21255472625327?)], 1)]

1+

1.00058518805631

var('s12,b1,b2,p') assume(b1>0) assume(b2>0) phi1(p)=(1/pi)*2^(3/2)*b1^(5/2)/(p^2+b1^2)^2 phi2(p)=(1/pi)*2^(3/2)*b2^(5/2)/(p^2+b2^2)^2 dv(p)=4*pi*p^2 s12=integral(phi1(p)*phi2(p)*dv(p),p,0,oo);s12

8*(pi*b1^3 - 3*pi*b1^2*b2 + 3*pi*b1*b2^2 -
pi*b2^3)*b1^(3/2)*b2^(3/2)/((b1^6 - 3*b1^4*b2^2 + 3*b1^2*b2^4 -
b2^6)*pi)

b1=1.001 b2=1.002 n( 8*(pi*b1^3 - 3*pi*b1^2*b2 + 3*pi*b1*b2^2 - pi*b2^3)*b1^(3/2)*b2^(3/2)/((b1^6 - 3*b1^4*b2^2 + 3*b1^2*b2^4 - b2^6)*pi) )

MomentumRepresentation

1 day ago by reibaretti