1/* Use the variational method to estimate the eigenvalues of
2
3   -f'' + (x^2 + epsilon * x^4) * f = mu * f
4
5for epsilon near zero.  The hamiltonian is
6*/
7
8ham (e)  := -diff(e, x, 2) + (x^2  + epsilon * x^4 ) * e;
9
10/*  Assume a trial solution psi that is a linear combination of n+1 even
11order Hermite polynomials times a Gaussian function. We'll need to
12assign a value to n, load orthopoly, and assign psi.
13*/
14
15n : 3;
16
17if get('orthopoly,'version) = 'false then load("orthopoly")$
18
19psi : sum(c[2*k] * hermite(2*k,x) * exp(-x^2 / 2),k,0,n) / %pi^(1/4)$
20
21/*  The denominator %pi^(1/4) makes the computation easier. Let
22vars be a list of the unknown c's.  Although the c's really aren't
23positive, we'll set assume_pos to true; doing so prevents
24Maxima from asking lots of questions about the signs of the
25c's.
26*/
27
28vars : makelist(c[ 2*i ],i,0,n)$
29
30assume_pos : true;
31
32/* Define the L2 inner product with the match fix operator
33<< , >>. Everything is  real, so we don't need a conjugate.
34*/
35
36matchfix("<<", ">>")$
37
38"<<" (f, g) := integrate(expand (f * g), x,-inf, inf)$
39
40/*  Minimize << psi, ham(psi) >> subject to the constraint
41<< psi, psi >> =1; let mu be the Lagrange multiplier.
42*/
43
44min_this : << psi, ham(psi) >> - mu * << psi, psi >>;
45
46eqs : makelist(diff(min_this,vars[ i ]),i,1,n+1)$
47
48/* The equations are linear and homogeneous in the c's.  Demand
49that the coefficient matrix is singular.
50*/
51m_det : determinant(coefmatrix(eqs, vars))$
52m_det : ratsimp(m_det)$
53
54/* Solve for mu as power series in epsilon.  Thus assume
55mu = cf[ 0] + cf[1] * epsilon + ... + cf[solve_ord] epsilon^solve_ord.
56*/
57
58solve_ord : 3;
59pows : makelist(epsilon^i,i,0,solve_ord)$
60unks : makelist(cf[ i ],i,0,solve_ord)$
61eq : ev(m_det, mu = unks . pows)$
62eq : taylor(eq, epsilon, 0, solve_ord)$
63eq : expand(eq)$
64eq : makelist(coeff(eq,epsilon,i),i,0,solve_ord)$
65ans : algsys(eq, unks)$
66
67for i : 1 thru length(ans) do (
68       ans[ i ] : map(rhs, ans[ i ]) . pows)$
69
70ans : reverse(ans);
71
72/* Look at the solution graphically.*/
73
74plot2d(ans, [epsilon,0,0.25]);
75
76/*  Let's solve the equations using allroots instead of the series method. */
77
78f(x,k) := part(sort(map('rhs, allroots(subst('epsilon=x,m_det)))),k);
79
80/* Compare the allroots solution to the series solution. */
81
82plot2d([ans[1], '(f(epsilon,1))], [epsilon,0.0,0.4]);
83plot2d([ans[2], '(f(epsilon,2))], [epsilon,0.0,0.4]);
84plot2d([ans[3], '(f(epsilon,3))], [epsilon,0.0,0.4]);
85plot2d([ans[4], '(f(epsilon,4))], [epsilon,0.0,0.4]);
86
87plot2d(['(f(epsilon,1)),'(f(epsilon,2)),'(f(epsilon,3)),'(f(epsilon,4))],[epsilon,0,0.4]);
88
89remfunction(ham,"<<",f);
90remvalue(n,psi,vars,min_this,eqs,m_det,solve_ord,pows,unks,eq,ans);
91assume_pos : false;
92
93
94/* Let's apply a variational method to the potential x^2 / 2 + x^4. We'll assume
95a trial wavefuction of the form qo * exp(-%alpha * abs(x)^(2*n) / 2) where the
96parameters are %alpha and n. See "Post-Gaussian variational method for quantum anharmonic
97oscillator," by Akihiro Ogura." */
98
99kill(all)$
100assume(qo > 0, %alpha > 0, n > 1/2)$
101f : qo * exp(-%alpha * abs(x)^(2*n) / 2);
1021 = integrate(f^2,x,minf,inf);
103solve(%,qo);
104f : subst(second(%), f);
105v : x^2 / 2 + x^4$
106ham(f) := -diff(f,x,2) / 2 + v * f$
107energy : integrate(f * ham(f),x,minf,inf);
108eqs : [diff(energy,n), diff(energy,%alpha)]$
109load(mnewton)$
110newtonepsilon : 1.0e-15$
111sol : mnewton(eqs,[n,%alpha],[1.1, 2.0]);
112subst(sol, energy);
113
114