3
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ Try Orthogonalize[your polys, Integrate[#1 #2, {x, -a, a}] &] $\endgroup$ – wuyudi Sep 9 at 5:18
  • $\begingroup$ 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. $\endgroup$ – yarchik Sep 9 at 6:14
6
$\begingroup$

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

| improve this answer | |
$\endgroup$
  • $\begingroup$ 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". $\endgroup$ – yawnoc Sep 9 at 7:17
  • 1
    $\begingroup$ Done. [Extra characters to satisfy the bots.] $\endgroup$ – Michael E2 Sep 9 at 13:27
  • $\begingroup$ @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]]; $\endgroup$ – DUO Sep 10 at 1:13

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.