# Equation to Matlab code

조회 수: 1(최근 30일)
Ajai Singh 2022년 8월 1일
댓글: Ajai Singh 2022년 8월 2일
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. ##### 댓글 수: 5표시 이전 댓글 수: 4숨기기 이전 댓글 수: 4
John D'Errico 2022년 8월 2일
편집: John D'Errico 2022년 8월 2일
@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.

댓글을 달려면 로그인하십시오.

### 채택된 답변

John D'Errico 2022년 8월 2일
편집: John D'Errico 2022년 8월 2일
(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:
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.
##### 댓글 수: 3표시 이전 댓글 수: 2숨기기 이전 댓글 수: 2
Ajai Singh 2022년 8월 2일
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.

댓글을 달려면 로그인하십시오.

### 범주

Find more on Linear Algebra in Help Center and File Exchange

R2021a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!