# Error in Creating Orthogonal Polynomials

I'm trying to create my own set of polynomials orthogonal to weight $$w(x)=x^{14}$$ on $$[-a,a]$$. My code:

N1=10;
Int0[n_,x_]=Integrate[x*(DUO[n])^2*x^14,x];
Int1[n_,x_]=Integrate[(DUO[n])^2*x^14,x];

DUO[0]=1;
DUO[1]=(x-Divide[Int0[0,N1]-Int0[0,-N1],Int1[0,N1]-Int1[0,-N1]])*DUO[0];

DUO[n_]:=DUO[n]=((x-Divide[Int0[n-1,N1]-Int0[n-1,-N1],Int1[n-1,N1]-Int1[n-1,-N1]])*DUO[n-1])-(Divide[Int1[n-1,N1]-Int1[n-1,-N1],Int1[n-2,N1]-Int1[n-2,-N1]]*DUO[n-2]);
DUO[2]


However, when evaluating DUO[2], I keep on getting random and non-sensical errors. I'm following the Gram-Schmidt Process, and I checked it over several times, and I can't find anything wrong in my implementation (other than it being horribly inefficient :) ). So could someone help me find what's wrong here?

• Try Orthogonalize[your polys, Integrate[#1 #2, {x, -a, a}] &] – wuyudi Sep 9 at 5:18
• Wellcome here. Yes, your implementation is inefficient. 1) Besides computing the integral of a square of a polynomial you also will need to integral the product of different ones. 2) You can directly compute the definite integral, not need to substitute the limits by hands. – yarchik Sep 9 at 6:14

The problem is that Int0 and Int1 are defined with = (Set) rather than := (SetDelayed), so the RHS evaluates immediately. This is problematic because DUO[n] does not yet explicitly depend on x. For example, your definition for Int0 becomes

Int0[n_, x_] = 1/16 x^16 DUO[n]^2


(to see this, run ?Int0).

Also we should directly evaluate the definite integral rather than take the difference between indefinite integrals. We define our inner product:

ip[f_, g_] := Integrate[f g x^14, {x, -a, a}];
ip[f_] := ip[f, f];


Then we implement (18), (19), and (24) in the linked mathworld page:

p[0] = 1;
p[1] = (x - ip[x p[0], p[0]] / ip[p[0]]) * p[0];
p[n_] := p[n] =
Subtract[
(x - ip[x p[n - 1], p[n - 1]] / ip[p[n - 1]]) * p[n - 1],
ip[p[n - 1]] / ip[p[n - 2]] * p[n - 2]
];


We get:

Table[{n, p[n]}, {n, 0, 4}] // TableForm


$$\begin{array}{cc} 0 & 1 \\ 1 & x \\ 2 & x^2-\frac{15 a^2}{17} \\ 3 & x \left(x^2-\frac{15 a^2}{17}\right)-\frac{4 a^2 x}{323} \\ 4 & x \left(x \left(x^2-\frac{15 a^2}{17}\right)-\frac{4 a^2 x}{323}\right)-\frac{289}{399} a^2 \left(x^2-\frac{15 a^2}{17}\right) \\ \end{array}$$

Check orthogonality:

Table[ip[p[m], p[n]], {m, 0, 4}, {n, 0, 4}] // TableForm


$$\begin{array}{ccccc} \frac{2 a^{15}}{15} & 0 & 0 & 0 & 0 \\ 0 & \frac{2 a^{17}}{17} & 0 & 0 & 0 \\ 0 & 0 & \frac{8 a^{19}}{5491} & 0 & 0 \\ 0 & 0 & 0 & \frac{8 a^{21}}{7581} & 0 \\ 0 & 0 & 0 & 0 & \frac{128 a^{23}}{3661623} \\ \end{array}$$

• Could someone with higher reputation please add the "check orthogonality" output? When I try to insert the LaTeX I keep getting "Your post appears to contain code that is not properly formatted as code". – yawnoc Sep 9 at 7:17
• Done. [Extra characters to satisfy the bots.] – Michael E2 Sep 9 at 13:27
• @user72028 Thank you for this this highly detailed answer--- +1 and accept. I also want to point out for future reference (since it took me like 3 hours to find this, and I don't want anyone going through the same effort), you can evaluate the indefinite integral at any point a like this: Int0[n_,a_]:=CompoundExpression[Expr:=Integrate[x*(DUO[n])^2*x^14,x],Func[x_]:=Evaluate[Expr],Func[a]]; – DUO Sep 10 at 1:13