Gaussian Quadrature ( Legendre Polynomials )

14 次查看（过去 30 天）
Ho Nam Ernest Yim
Ho Nam Ernest Yim 2017-11-13
评论： Karan Gill 2017-11-13
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

请先登录，再回答此问题。

采纳的回答

John D'Errico
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!

Translated by