share/orthopoly/variational_method.dem

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