481 - Legendre Polynomials and Gaussian Quadrature

738 days ago by Professor481

# Let's solve the 2-pt GQ parameters "by hand" var('x1,x2,w1,w2') solve([w1+w2-2, w1*x1+w2*x2, w1*x1^2+w2*x2^2 - 2/3, w1*x1^3+w2*x2^3],(x1,x2,w1,w2)) 
       
[[x1 == 1/3*sqrt(3), x2 == -1/3*sqrt(3), w1 == 1, w2 == 1], [x1 ==
-1/3*sqrt(3), x2 == 1/3*sqrt(3), w1 == 1, w2 == 1]]
[[x1 == 1/3*sqrt(3), x2 == -1/3*sqrt(3), w1 == 1, w2 == 1], [x1 == -1/3*sqrt(3), x2 == 1/3*sqrt(3), w1 == 1, w2 == 1]]
# Technically, this should find the 3-pt GQ parameters but Sage can't do it var('x1,x2,x3,w1,w2,w3') solve([w1+w2+w3-2, w1*x1+w2*x2+w3*x3, w1*x1^2+w2*x2^2+w3*x3^2 - 2/3, w1*x1^3+w2*x2^3+w3*x3^3, w1*x1^4+w2*x2^4+w3*x3^4 - 2/5, w1*x1^5+w2*x2^5+w3*x3^5],(x1,x2,x3,w1,w2,w3)) 
       
^C
[w1 + w2 + w3 - 2, w1*x1 + w2*x2 + w3*x3, w1*x1^2 + w2*x2^2 + w3*x3^2 -
2/3, w1*x1^3 + w2*x2^3 + w3*x3^3, w1*x1^4 + w2*x2^4 + w3*x3^4 - 2/5,
w1*x1^5 + w2*x2^5 + w3*x3^5]
__SAGE__
^C
[w1 + w2 + w3 - 2, w1*x1 + w2*x2 + w3*x3, w1*x1^2 + w2*x2^2 + w3*x3^2 - 2/3, w1*x1^3 + w2*x2^3 + w3*x3^3, w1*x1^4 + w2*x2^4 + w3*x3^4 - 2/5, w1*x1^5 + w2*x2^5 + w3*x3^5]
__SAGE__
deg = 5 print("The roots of Legendre polynomial of degree "+str(deg)) solns = solve(legendre_P(deg,x),x,solution_dict=True) print(solns) print("are the roots of P(x) = "+str(legendre_P(deg,x))) 
       
The roots of Legendre polynomial of degree 5
[{x: -1/3*sqrt(2/7*sqrt(70) + 5)}, {x: 1/3*sqrt(2/7*sqrt(70) + 5)}, {x:
-1/3*sqrt(-2/7*sqrt(70) + 5)}, {x: 1/3*sqrt(-2/7*sqrt(70) + 5)}, {x: 0}]
are the roots of P(x) = 63/8*x^5 - 35/4*x^3 + 15/8*x
The roots of Legendre polynomial of degree 5
[{x: -1/3*sqrt(2/7*sqrt(70) + 5)}, {x: 1/3*sqrt(2/7*sqrt(70) + 5)}, {x: -1/3*sqrt(-2/7*sqrt(70) + 5)}, {x: 1/3*sqrt(-2/7*sqrt(70) + 5)}, {x: 0}]
are the roots of P(x) = 63/8*x^5 - 35/4*x^3 + 15/8*x
xx = [] for k in range(deg): xx.append(solns[k][x]) show(xx[k]) 
       

                                
                            

                                
# Next let's get weights Pprime = derivative(legendre_P(deg,x),x) print("The "+str(deg)+"th Legendre derivative:") Pprime.show() print("yields x-values and quadrature coefficients:") w = [] for k in range(deg): w.append(2/((1-xx[k]^2)*Pprime(x = xx[k])^2)) show([xx[k],w[k]]) 
       
The 5th Legendre derivative:
yields x-values and quadrature coefficients:

                                
                            
The 5th Legendre derivative:
yields x-values and quadrature coefficients:

                                
# so the associated gaussian quadrature rule that integrates on [-1,1] exactly a polynomial of highest degree possible f = x^6 Integral = 0 for k in range(deg): print(xx[k]) Integral += w[k]*f(x=xx[k]) show(Integral.n()) show(integrate(f,x,-1,1)) 
       
-1/3*sqrt(2/7*sqrt(70) + 5)
1/3*sqrt(2/7*sqrt(70) + 5)
-1/3*sqrt(-2/7*sqrt(70) + 5)
1/3*sqrt(-2/7*sqrt(70) + 5)
0

                                
                            
-1/3*sqrt(2/7*sqrt(70) + 5)
1/3*sqrt(2/7*sqrt(70) + 5)
-1/3*sqrt(-2/7*sqrt(70) + 5)
1/3*sqrt(-2/7*sqrt(70) + 5)
0

                                
# using the 2-pt gaussian quadrature rule (5/2*(sin((-5/sqrt(3)+11)/2)+sin((5/sqrt(3)+11)/2))).n() 
       
-0.448286680695642
-0.448286680695642
# General conversion from [-1,1] to [a,b] # u = -1*(x-b)/(a-b) + 1*(x-a)/(b-a) = (2x - (a+b))/(b-a) # x = ( (b-a)*u + (a+b) )/2 # dx = (b-a)/2 du 
       
f = sin(x) a = 3 b = 4 GQ = (b-a)/2 * ( f( x = (( (b-a)*xx[0] + (a+b) )/2)) + f( x = (((b-a)*xx[1] + (a+b) )/2)) ) GQ.n() 
       
-0.324801704989511
-0.324801704989511
f = sin(x) a = 3 b = 6 # Below is what would be in the function that you would call to do Gaussian Quadrature on pieces # of a larger integral GQ = 0 for k in range(deg): GQ += w[k]*f( x = (( (b-a)*xx[k] + (a+b) )/2)) GQ = GQ*(b-a)/2 pretty_print(GQ.n()) integrate(f,(x,a,b)).n() 
       
-1.95016278325081
-1.95016278325081