Gaussian Quadrature ( Legendre Polynomials )
조회 수: 15 (최근 30일)
이전 댓글 표시
I have tried to create a function that will find the coefficients of the nth order Legendre polynomial without using syms x, but I have got stuck on how to actually get the coefficients with 2 unknowns in my equation... Can anyone help a bit??
Here's my code :
function [y] = LegendrePoly(n)
if n == 0
y = 1;
elseif n == 1
y = x;
else
x = zeros(n+1,1);
x(n+1) = 1;
x_1 = zeros(n+1,1);
x_1(n) = 1;
for i=2:n;
y(i)=((2*x_1(n))/x_1(n))*x*y(i-1)-(x(n-1)/x_1(n))*y(i-2);
end
end
댓글 수: 1
Karan Gill
2017년 11월 13일
There's an example on this topic: https://www.mathworks.com/help/symbolic/examples/gauss-laguerre-quadrature-evaluation-points-and-weights.html. Is it helpful?
채택된 답변
John D'Errico
2017년 11월 13일
편집: John D'Errico
2017년 11월 13일
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?
y = x;
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); % recurrence
[P_nm1,P_n] = deal(P_n,P_np1); % shift
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
coefs =
1.5 0 -0.5
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
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Polynomials에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!