1618.sagews

# MATHEMATICAL METHODS OF Science and Engineering
# Lectures on Quantum Mechanics /Density Functional theory
#http://www.uprh.edu/rbaretti
#http://www1.uprh.edu/rbaretti/Methodsoftheoreticalphysics.htm
# http://www1.uprh.edu/rbaretti/methodsoftheoreticalphysicspart2.htm
# http://www1.uprh.edu/rbaretti/methodsoftheoreticalphysicspart3.htm
# http://www1.uprh.edu/rbaretti/methodsoftheoreticalphysicspart4.htm
# http://www1.uprh.edu/rbaretti/methodsoftheoreticalphysicsPart5.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 (most recent call last): norm=integral_numerical(phi1s(p)*phi2s(p)*dv(p),0,1.5,max-points=500)[0] File "", line 1, in <module> File "/tmp/tmpVI6ZIf/___code___.py", line 8 norm=integral_numerical(phi1s(p)*phi2s(p)*dv(p),_sage_const_0 ,_sage_const_1p5 ,max-points=_sage_const_500 )[_sage_const_0 ] 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

1

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 (most recent call last): File "<stdin>", line 1, in <module> File "_sage_input_5.py", line 9, in <module> exec compile(ur'open("___code___.py","w").write("# -*- coding: utf-8 -*-\n" + _support_.preparse_worksheet_cell(base64.b64decode("QT1tYXRyaXgoW1sgLTAuODUzNDhFKzAzICwtMC4xNTA2MEUrMDIgIF0sWyAtMC4xNTA2MEUrMDIgLCAgLTAuMjMxMjlFKzAzIF1dKQpBLmVpZ2VudmVjdG9ycygp"),globals())+"\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single') File "", line 1, in <module> File "/tmp/tmpsooe2n/___code___.py", line 4, in <module> exec compile(ur'A.eigenvectors()' + '\n', '', 'single') File "", line 1, in <module> File "element.pyx", line 305, in sage.structure.element.Element.__getattr__ (sage/structure/element.c:2743) File "parent.pyx", line 237, in sage.structure.parent.getattr_from_other_class (sage/structure/parent.c:2844) File "parent.pyx", line 169, in sage.structure.parent.raise_attribute_error (sage/structure/parent.c:2611) 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 (most recent call last): sage: A.eigenvalues() File "", line 1, in <module> File "/tmp/tmppubkUF/___code___.py", line 3 http://sagemath.org/doc/constructions/linear_algebra.html ^ 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)   )