481 - Legendre Polynomials and Gaussian Quadrature -- Sage

# 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])

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{3} \, \sqrt{\frac{2}{7} \, \sqrt{70} + 5}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{3} \, \sqrt{\frac{2}{7} \, \sqrt{70} + 5}$

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{3} \, \sqrt{-\frac{2}{7} \, \sqrt{70} + 5}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{3} \, \sqrt{-\frac{2}{7} \, \sqrt{70} + 5}$

$\newcommand{\Bold}[1]{\mathbf{#1}}0$

# 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:

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{315}{8} \, x^{4} - \frac{105}{4} \, x^{2} + \frac{15}{8}$

yields x-values and quadrature coefficients:

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{1}{3} \, \sqrt{\frac{2}{7} \, \sqrt{70} + 5}, -\frac{16003008}{25 \, {\left({\left(2 \, \sqrt{70} + 35\right)}^{2} - 84 \, \sqrt{70} - 1281\right)}^{2} {\left(\sqrt{70} - 14\right)}}\right]$

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\frac{1}{3} \, \sqrt{\frac{2}{7} \, \sqrt{70} + 5}, -\frac{16003008}{25 \, {\left({\left(2 \, \sqrt{70} + 35\right)}^{2} - 84 \, \sqrt{70} - 1281\right)}^{2} {\left(\sqrt{70} - 14\right)}}\right]$

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[-\frac{1}{3} \, \sqrt{-\frac{2}{7} \, \sqrt{70} + 5}, \frac{16003008}{25 \, {\left({\left(2 \, \sqrt{70} - 35\right)}^{2} + 84 \, \sqrt{70} - 1281\right)}^{2} {\left(\sqrt{70} + 14\right)}}\right]$

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\frac{1}{3} \, \sqrt{-\frac{2}{7} \, \sqrt{70} + 5}, \frac{16003008}{25 \, {\left({\left(2 \, \sqrt{70} - 35\right)}^{2} + 84 \, \sqrt{70} - 1281\right)}^{2} {\left(\sqrt{70} + 14\right)}}\right]$

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, \frac{128}{225}\right]$

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

$\newcommand{\Bold}[1]{\mathbf{#1}}0.285714285714285$

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{2}{7}$

-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()

$\newcommand{\Bold}[1]{\mathbf{#1}}-1.95016284845677$

-1.95016278325081

-1.95016278325081

481 - Legendre Polynomials and Gaussian Quadrature

749 days ago by Professor481