Error in Creating Orthogonal Polynomials

Question

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?

yawnoc · Accepted Answer

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

Error in Creating Orthogonal Polynomials

One Answer

Add your own answers!

Ask a Question