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