2
$\begingroup$

Regarding to $e^{-ax-bx^2}$ as the weight function , how can I generate a set of orthogonal polynomials? I just know the Orthogonalize command in Mathematica which seems that uses 1 as weight function by Gram-Schmidt method, other commands such as HermiteH and LaguerreL generate special polynomials. But I don't know how to do my case. In fact I try to perform the following step which is a part of one article:

Generate a set of orthogonal polynomials with the square of the Ground-state wave function as the weight function (WF) adopting any of the standard algorithms

where Ground-state wave function is $e^{-ax-bx^2}$.

$\endgroup$
  • $\begingroup$ What inner product are you using for these polynomials? $\endgroup$ – Marius Ladegård Meyer Jul 4 at 12:58
1
$\begingroup$

If you remove the $x^n$ in your $A$ it will work. The ground state wave-function $\psi(x)=\exp(-ax-bx^2)$ and the weight function is its square $\psi^2(x)=\exp(-2 a x-2 b x^2)$. It helps avoid numerical instability later on if we normalize the weight function first by doing: $$ \hat{w}(x)=\psi^2(x)*\left(\int_{-\infty}^\infty \psi^2(x)\mathrm{d}x\right)^{-1} =\frac{\psi^2(x)\sqrt{b}}{\sqrt{\frac{\pi }{2}} e^{\frac{a^2}{2 b}}} $$

(* ground state wave function *)
ψ = Exp[-a x - b x^2];

(* weight function is the square of ψ *)
w = ψ^2;

(* norm of the weight function *)
normw = Integrate[w, {x, -∞, ∞}];

(* normalized weight function *)
wnormed = w/normw;

(* orthogonalizes the monomials 1,x,x^2,x^3,...,x^n assuming an inner-product convolution with normalized weight fn *)
symbasis[n_?IntegerQ] :=
 Orthogonalize[x^# & /@ Range[0, n], 
  Integrate[wnormed * #1 * #2, {x, -∞, ∞}] &]

For example, you can then get a basis (in symbolic $a,b$) for $n=3$ with symbasis[3]. Replace $a,b$ with numerical values if necessary e.g: symbasis[3] /. {a->0.99985, b->0.00015}. This process of getting a symbolic basis then replacing the values is slower than numerical integration but it seems to avoid the problems with precision. With this replacement we get a basis:

basis = symbasis[3] /. {a->0.99985, b->0.00015}
(* results: 
{1,
 0.0244949 (3332.83 + x), 
 0.000424264 (-1.11094*10^7 + x^2 + 6665.67 (3332.83 + x)), 
 6.*10^-6 (3.7037*10^10 + x^3 - 3.33283*10^7 (3332.83 + x) + 9998.5 (-1.11094*10^7 + x^2 + 6665.67 (3332.83 + x)))
}*)

We can verify the basis functions are orthogonal as follows:

(* all six pairs of the 4 basis functions are orthogonal *)
Integrate[#[[1]]*#[[2]]*(wnormed /. {a -> 0.99985, b -> 0.00015}) , {x, -∞, ∞}] & /@ Subsets[basis, {2}]

(* result: {0., 0., 0., 0., 0., 0.} *)
| improve this answer | |
$\endgroup$
  • $\begingroup$ Dear @flinty thanks a lot for your reply. I ran it in Mathematica and I got a list with 3 elements. Namely it generated a set of orthogonal polynomials (with 3 members) just using above function as weight function? However it's hard for me to understand your code. Can you explain details of your code? Also I added the original sentences appeared in the paper. $\endgroup$ – Wisdom Jul 4 at 14:02
  • $\begingroup$ Yes that's exactly what it does. I'm afraid I don't do much physics though. $\endgroup$ – flinty Jul 4 at 14:06
  • $\begingroup$ Dear @flinty Good, no matter because this is a purely mathematical step. However I want to understand the details of your code, I'm new in Mathematica. Can you tell me every line of your code what does? Thanks again. $\endgroup$ – Wisdom Jul 4 at 14:09
  • $\begingroup$ You read the quote which I included in my question? @flinty $\endgroup$ – Wisdom Jul 4 at 14:12
  • $\begingroup$ I have commented the code. What specifically do you not understand? $\endgroup$ – flinty Jul 4 at 14:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.