# 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}] &] Sep 9, 2020 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. Sep 9, 2020 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". Sep 9, 2020 at 7:17
• Done. [Extra characters to satisfy the bots.] Sep 9, 2020 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]]; Sep 10, 2020 at 1:13