Equation to Matlab code

Can someone please help to me understand how the equation (eq. 47) and code given in the picture relate ?

Any hint or suggestion would be of great help.

In the equation c is the order of the polynomial.

Thank you.

[Rik]

There are several algebraic steps between the two forms. The x seems to be substituted already.

I must admit I don't understand the mathematical description well enough to understand why the code is equivalent.

[Torsten]

I guess LPC(i,j) is the value of j-th coefficient of the i-th Legendre polynomial.

[jessupj]

i think OP needs to define what the intended domain is. there's no x anywhere here.

start with

x=linspace(0,1,25);

[Jeffrey Clark]

@Ajai Singh, @Torsten is correct - the LPC values are the computed multipliers for the x^(c-n) terms of the Legendre polynomials. So LPC(1,1) is the multiplier for the zero order Legendre polynomial, which is 1*x^0; similarly LPC(2,1:2) are the mulpliers for the 1st order polynomial 0*x^0 + 1*x^1. Check that the MATLAB code provided does compute the 2nd order polynomial multipliers when c = 3 to be LPC(3,1:3) = [-0.5 0 1.5] as in -0.5*x^0 + 0*x^1 + 1.5*x^2.

A table for the first 11 (0->10) can be seen in this article Legendre polynomials - Wikipedia

Note that MATLAB functions that operate on this short hand description of polynomials generally expects them in reverse order. See Polynomial roots - MATLAB roots (mathworks.com)

[John D'Errico]

@jessupj - The Legendre polynomials from this family will be orthognal on the interval [-1,1], NOT [0,1]. In fact though, they could be evaluated on ANY interval, as they are still just polynomials. Some people choose to define the Legendre polynomials on the interval [0,1], which is entirely doable. But were you to do that, you would find a different family of polynomials. Regardless, you do not need to define what x is before you create any polynomials.

Answer 1

[John D'Errico]

(I really hope I am not doing your homework here, but the question is not posed as if it is, even though I assume you are a student.)

My take? Calling that a translation is a fairly liberal use of the word. I would also point out that it orders the coefficients of the polynomials in the wrong way. The MATLAB standard is to start with the highest order term FIRST.

On the interval [-1,1] the family of Legendre polynomials are an orthogonal family. We can find them listed in this link:

https://en.wikipedia.org/wiki/Legendre_polynomials

You can see the same 3 term recurrence relation posed in your question as Bonnet's recursion formula in that link.

The first few such polynomials are:

p_0(x) = 1

p_1(x) = x

p_2(x) = (3*x^2 - 1)/2

As you can see, the first two rows of LPC give us the first two polynomials, since we will interpret the coefficients in those rows as the coefficients of the powers of x. Essentially, those first two rows start the recursion.

But now, suppose I decided to write a truly direct translation of the recurrence relation? What would it look like? I might do this:

C = 7;
LPC = zeros(C+1,C+1);
LPC(1,1) = 1;
LPC(2,1:2) = [0 1];
for c = 1:C-1
  LPC(c+2,1:c+2) = ((2*c+1)*[0,LPC(c+1,1:c+1)] - c*[LPC(c,1:c),0 0])/(c+1);
end

For example, do you see that multiplying the polynomial represented by coefficients [0 1], by x is represented by the code fragment where we shift the coefficients over by pre-pending a zero to the vector? So [0 1] to represent x would turn into [0 0 1], representing x^2.

LPC

LPC = 8×8

1.0000 0 0 0 0 0 0 0 0 1.0000 0 0 0 0 0 0 -0.5000 0 1.5000 0 0 0 0 0 0 -1.5000 0 2.5000 0 0 0 0 0.3750 0 -3.7500 0 4.3750 0 0 0 0 1.8750 0 -8.7500 0 7.8750 0 0 -0.3125 0 6.5625 0 -19.6875 0 14.4375 0 0 -2.1875 0 19.6875 0 -43.3125 0 26.8125

We can make that prettier and more easily read using syms.

syms x
LPC*x.^(0:C).'

ans = 

As you can see, this does generate the table of Legendre polynomials. We can see the same table in the wiki link I gave. And the code I wrote is a fairly faithful translation of the recurrence relation. It should be pretty easy to read, once you understand that multiplying a polynomial by x is equivalent to shifting the coefficients over and pre-pending a zero.

Perhaps a more interesting question is, why does the code you show also create the same set of coefficients? An analogy here would be, if I gave you the formula for a convolution of two functions, as a convolution. But suppose I wrote code using a Fourier transform, where I take the transforms of each function, multiply them in the Fourier domain, and then take the inverse Fourier transform? The result would still be the convolution of the functions, but I would surely not call it a literal translation of a convolution.

Anyway, why does the code written in your question work? You should see that I used only a single loop. Inside the loop, I extracted the coefficients of an entire lower order polynomial at once, then used them directly as a vector of coefficients. The loop you had is a double loop, but I'll claim it is a terribly poor implementation, even at that. Inside the inner loop, they first define LPC(i,j+1), and THEN define LPC(i,j), each time in the loop. UGH. Sorry, but that is just a convoluted code. What it does is allow you to not need a special case for the constant coefficient. I'll re-write my own code, now using a double, nested loop, DIRECTLY from the formula. As you should see, I create the constant term outside the loop.

C = 7;
LPC = zeros(C+1,C+1);
LPC(1,1) = 1;
LPC(2,1:2) = [0 1];
for c = 1:C-1
  LPC(c+2,1) = -c*LPC(c,1)/(c+1); % the constant coefficient is assigned outside the loop
  for j = 1:c+1
    LPC(c+2,j+1) = ((2*c+1)*LPC(c+1,j) - c*LPC(c,j+1))/(c+1);
  end
end

LPC*x.^(0:C).'

ans = 

As you should see, this does generate the same set of polynomials, using a nested loop.

[Ajai Singh]

Thank you @John D'Errico . I figured it out. No it was not a home work problem. Part of my research was stuck here as I am was trying to relate the mathematical equation to its matlab code. Turns out we got bit lost but thanks for clearing the doubt.

[John D'Errico]

The code you had was difficult to read, because while it was working, it did the job in a strange way. The simpler approach, where I deal with the complete polynomials in a single loop as vectors of polynomial coefficients is I think easier to follow. Easy to read code is always the best code, because then you can read it and debug it and most importantly, understand it.

[Ajai Singh]

I agree , we referred to a papers supplimentary material and got that piece of code. The way i figured it out was comparing the code and the given equation and then solved c in terms of i and then everything was clear. Thank you for the suggestion. It made the code more readable. I appreciate it.