NO NO NO. You say that you don't want to use syms here. So why in the name of god and little green apples are you doing things like this?
All you need to store are the COEFFICIENTS!!!!!!! It is a polynomial. You have not passed x into the function to evaluate. That you can do trivially with polyval anyway.
Set up LegendrePoly to return a vector of length n+1. The coefficients of that polynomial.
I'm feeling too lazy to crack open Abramowitz & Stegun, so this page gives the recurrence relation and the first two polynomials.
https://en.wikipedia.org/wiki/Legendre_polynomials
Thus, P_0 = 1. P_1 = x, and:
(n+1)*P_(n+1) = (2*n+1)*P_n*x - n*P_(n-1)
Next, DON'T EVER write recurrence relations like this recursively, unless you use a memoized form. The problem is that the number of recursive calls grows exponentially. It is far simpler to just use a loop anyway.
function coefs = LegendrePoly(n)
if n == 0
coefs = 1;
elseif n == 1
coefs = [1 0];
else
P_nm1 = 1;
P_n = [1 0];
for i=1:(n-1);
P_np1 = ((2*i+1)*[P_n,0] - i*[0,0,P_nm1])/(i+1);
[P_nm1,P_n] = deal(P_n,P_np1);
end
coefs = P_np1;
end
Think about what [P_n,0] does, in terms of polynomial coefficients.
Thus we know that P_2=(3*x^2-1)/2. For n==2, we get
Likewise, for n==5, we get
coefs
coefs =
7.875 0 -8.75 0 1.875 0
From the table of polynomials, I got what I expected.
[63/8 0 -70/8 0 15/8 0]
ans =
7.875 0 -8.75 0 1.875 0
You can evaluate that 5th degree polynomial simply as:
polyval(coefs,0.5)
ans =
0.08984375
Finally, in order to use them as polynomials for Gaussian quadrature, you will need the derivative polynomials too.
polyder(coefs)
ans =
39.375 0 -26.25 0 1.875
Its been a while since I had to derive the Gaussian quadrature but you need some roots too. Again, trivial.
roots(coefs)
ans =
0
-0.906179845938664
-0.538469310105683
0.906179845938664
0.538469310105683
Best Answer