# How to generate a set of orthogonal polynomials with a special weight function

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}$$.

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

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. Replace $$a,b$$ with numerical values if necessary e.g: symbasis /. {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 /. {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[#[]*#[]*(wnormed /. {a -> 0.99985, b -> 0.00015}) , {x, -∞, ∞}] & /@ Subsets[basis, {2}]

(* result: {0., 0., 0., 0., 0., 0.} *)

• 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. – Wisdom Jul 4 at 14:02
• Yes that's exactly what it does. I'm afraid I don't do much physics though. – flinty Jul 4 at 14:06
• 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. – Wisdom Jul 4 at 14:09
• You read the quote which I included in my question? @flinty – Wisdom Jul 4 at 14:12
• I have commented the code. What specifically do you not understand? – flinty Jul 4 at 14:34