381 - A11 - Numerical Integration -- Sage

John Travis

Mississippi College

Random lecture demonstrations dealing with numerical Integration.

Starting with Newton-Coates formulas...

For numerical integration, one can create approximations to integrals (and derivatives) by simply taking the integrand and approximating it with a polynomial (such as Lagrange) and then integrating the interpolating polynomial exactly...since it is only just a polynomial.

# Approximate an ugly function with a quadratic and compare integrals and derivative p = 79/56*x^2-433/56*x+303/28 f = 9/(1+x^3) G = plot(p,x,1,3) G += plot(f,x,1,3,color='red',title='Integration') fp = derivative(f,x) pp = derivative(p,x) H = plot(fp,x,1,3, color='red',title='Differentiation') H += plot(pp,x,1,3) show(G) show(H) I = integrate(9/(1+x^3),x,1,3) pretty_print(n(I)) N = integrate(79/56*x^2-433/56*x+303/28,x,1,3) pretty_print(n(N))

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

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

Consider the Composite Trapezoidal Rule

# Using an example f = sin(x) a = 0 b = 3/4*pi n = 32 delx = (b-a)/n xs = [] for k in range(n+1): xs.append(a+k*delx) end G = plot(f,x,a,b) for k in range(n+1): G += point((xs[k],f(x=xs[k])),color='red') end show(G) exact = integrate(f,x,a,b) trap = (f(x=xs[0])+f(x=xs[n]))/2 for k in range(1,n): trap += f(x=xs[k]) end trap = trap*delx err = exact-trap pretty_print(exact.n()) pretty_print(trap.n()) pretty_print(err.n())

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

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

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

For Simpson's rule, approximate the actual integrand with a quadratic evaluated at equally spaced values from a to b and integrate that quadratic exactly. Notice, this will only require that the integrand be evaluated at three pionts over the interval.

reset() var('a,b,h,x0,x1,x2,y0,y1,y2,h') x0 = a x1 = (a+b)/2 x2 = b assume(a<b) h=(b-a)/2 L0 = (x-x1)*(x-x2)/((x0-x1)*(x0-x2)) show(L0) L1 = (x-x0)*(x-x2)/((x1-x0)*(x1-x2)) show(L1) L2 = (x-x0)*(x-x1)/((x2-x0)*(x2-x1)) show(L2) I0 = integrate(L0,x,a,b).factor() I1 = integrate(L1,x,a,b).factor() I2 = integrate(L2,x,a,b).factor() show(I0) show(I1) show(I2) I = I0*y0+I1*y1+I2*y2 show(I)

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(a + b - 2 \, x\right)} {\left(b - x\right)}}{{\left(a - b\right)}^{2}}$

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{4 \, {\left(a - x\right)} {\left(b - x\right)}}{{\left(a - b\right)}^{2}}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(a + b - 2 \, x\right)} {\left(a - x\right)}}{{\left(a - b\right)}^{2}}$

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{6} \, a + \frac{1}{6} \, b$

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{2}{3} \, a + \frac{2}{3} \, b$

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{6} \, a + \frac{1}{6} \, b$

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{6} \, {\left(a - b\right)} y_{0} - \frac{2}{3} \, {\left(a - b\right)} y_{1} - \frac{1}{6} \, {\left(a - b\right)} y_{2}$

So, this suggests "Simpson's Rule"

$\int_a^b f(x) dx \approx (b-a)\left [ \frac{1}{6}f(a) + \frac{4}{6}f(\frac{a+b}{2}) + \frac{1}{6}f(b) \right ]$

# Using an example f = sin(x) var('a,b,h,x0,x1,x2,y0,y1,y2,h') x0 = a x1 = (a+b)/2 x2 = b aa = 0 bb = 3*pi/4 G = plot(f,x,aa,bb) y0 = f(x=aa) y1 = f(x=(aa+bb)/2) y2 = f(x=bb) L0 = (x-x1)*(x-x2)/((x0-x1)*(x0-x2)) L1 = (x-x0)*(x-x2)/((x1-x0)*(x1-x2)) L2 = (x-x0)*(x-x1)/((x2-x0)*(x2-x1)) p = L0(a=aa,b=bb)*y0 + L1(a=aa,b=bb)*y1 + L2(a=aa,b=bb)*y2 G += plot(p,x,aa,bb,color='red') show(G) pretty_print(html("$\int f(x) dx =$"+str(integrate(f,x,aa,bb).n()))) pretty_print(html("$\int p(x) dx =$"+str(integrate(p,x,aa,bb).n()))) delx = (bb-aa) Simp = (delx*((1/6)*y0 + (4/6)*y1 + (1/6)*y2)).n() pretty_print('And just using the three y-values with coefficients from above: ',Simp)

$\int f(x) dx =$ 1.70710678118655

$\int p(x) dx =$ 1.72890675970461

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|And|\phantom{\verb!x!}\verb|just|\phantom{\verb!x!}\verb|using|\phantom{\verb!x!}\verb|the|\phantom{\verb!x!}\verb|three|\phantom{\verb!x!}\verb|y-values|\phantom{\verb!x!}\verb|with|\phantom{\verb!x!}\verb|coefficients|\phantom{\verb!x!}\verb|from|\phantom{\verb!x!}\verb|above:| 1.72890675970461$

1.70710678118655

1.72890675970461

What about using a cubic to approximate.....giving "Simpson's 3/8 Rule"

The following needs to be edited...

$\int_a^b f(x) dx \approx (b-a)\left [ \frac{1}{6}f(a) + \frac{4}{6}f(\frac{a+b}{2}) + \frac{1}{6}f(b) \right ]$

var('a,b,h,x0,x1,x2,x3,y0,y1,y2,y3,h') assume(a<b) h=(b-a)/3 x0 = a x1 = a + h x2 = a + 2*h x3 = b L0 = (x-x1)*(x-x2)*(x-x3)/((x0-x1)*(x0-x2)*(x0-x3)) show(L0) L1 = (x-x0)*(x-x2)*(x-x3)/((x1-x0)*(x1-x2)*(x1-x3)) show(L1) L2 = (x-x0)*(x-x1)*(x-x3)/((x2-x0)*(x2-x1)*(x2-x3)) show(L2) L3 = (x-x0)*(x-x1)*(x-x2)/((x3-x0)*(x3-x1)*(x3-x2)) show(L3) I0 = integrate(L0,x,a,b).factor() I1 = integrate(L1,x,a,b).factor() I2 = integrate(L2,x,a,b).factor() I3 = integrate(L3,x,a,b).factor() show(I0) show(I1) show(I2) show(I3) I = I0*y0+I1*y1+I2*y2+I3*y3 show(I)

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{{\left(2 \, a + b - 3 \, x\right)} {\left(a + 2 \, b - 3 \, x\right)} {\left(b - x\right)}}{2 \, {\left(a - b\right)}^{3}}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{9 \, {\left(a + 2 \, b - 3 \, x\right)} {\left(a - x\right)} {\left(b - x\right)}}{2 \, {\left(a - b\right)}^{3}}$

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{9 \, {\left(2 \, a + b - 3 \, x\right)} {\left(a - x\right)} {\left(b - x\right)}}{2 \, {\left(a - b\right)}^{3}}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{{\left(2 \, a + b - 3 \, x\right)} {\left(a + 2 \, b - 3 \, x\right)} {\left(a - x\right)}}{2 \, {\left(a - b\right)}^{3}}$

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{8} \, a + \frac{1}{8} \, b$

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{3}{8} \, a + \frac{3}{8} \, b$

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{8} \, a + \frac{1}{8} \, b$

$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{1}{8} \, {\left(a - b\right)} y_{0} - \frac{3}{8} \, {\left(a - b\right)} y_{1} - \frac{3}{8} \, {\left(a - b\right)} y_{2} - \frac{1}{8} \, {\left(a - b\right)} y_{3}$

This suggests "Simpson's 3/8 Rule"

$\int_a^b f(x) dx \approx (b-a)\left [ \frac{1}{8}f(a) + \frac{3}{8}f(a + \frac{b-a}{3}) + \frac{3}{8}f(a + 2\frac{b-a}{3}) + \frac{1}{8}f(b) \right ]$

# Using an example f = sin(x) var('a,b,h,x0,x1,x2,x3,y0,y1,y2,y3,h') assume(a<b) h=(b-a)/3 x0 = a x1 = a + h x2 = a + 2*h x3 = b aa = 0 bb = 3*pi/4 G = plot(f,x,aa,bb) y0 = f(x=aa) y1 = f(x=aa + (bb-aa)/3) y2 = f(x=aa + 2*(bb-aa)/3) y3 = f(x=bb) L0 = (x-x1)*(x-x2)*(x-x3)/((x0-x1)*(x0-x2)*(x0-x3)) L1 = (x-x0)*(x-x2)*(x-x3)/((x1-x0)*(x1-x2)*(x1-x3)) L2 = (x-x0)*(x-x1)*(x-x3)/((x2-x0)*(x2-x1)*(x2-x3)) L3 = (x-x0)*(x-x1)*(x-x2)/((x3-x0)*(x3-x1)*(x3-x2)) p = L0(a=aa,b=bb)*y0 + L1(a=aa,b=bb)*y1 + L2(a=aa,b=bb)*y2 + L3(a=aa,b=bb)*y3 G += plot(p,x,aa,bb,color='red') show(G) pretty_print(html("$\int f(x) dx =$"+str(integrate(f,x,aa,bb).n()))) pretty_print(html("$\int p(x) dx =$"+str(integrate(p,x,aa,bb).n()))) delx = (bb-aa) Simp = delx*((1/8)*y0 + (3/8)*y1 + (3/8)*y2 + (1/8)*y3) pretty_print('And just using the three y-values with coefficients from above: ',Simp.n())

$\int f(x) dx =$ 1.70710678118655

$\int p(x) dx =$ 1.71661348472682

1.70710678118655

1.71661348472682

Since each of these formulas are scaled by the interval width b-a and only use a few points on that interval then if the interval is long then the approximations might be less accurate.  Breaking up a larger interval therefore into pieces and reusing each formula over and over on the smaller pieces might be advantageous.  Adding up the bits give the final result.

# Simpson's Composite Rule var('t') # f = sin(x) f = 1/(1+x^3) aa = 0 bb = 3*pi/4 G = plot(f,x,aa,bb) n = 24 # break the entire integral up into n pieces and do each one with the rule above and then accumulate. Make this even so that we can approximate with n/2 quadratics. delx = (bb-aa)/n xs = [] for k in range(n+1): temp = aa+k*delx xs.append(temp) G += points((temp,f(x=temp)),color='red') end ys = [] for k in range(n+1): ys.append(f(x=xs[k]).n()) end S = ys[0] + ys[n] G += text(1,(xs[0],f(x=xs[0])+0.1),color='green') G += text(1,(xs[n],f(x=xs[n])+0.1),color='green') for k in range(1,n/2+1): temp = 2*k-1 S += 4*ys[temp] G += text(4,(xs[temp],f(x=xs[temp])+0.1),color='green') end for k in range(1,n/2): temp = 2*k S += 2*ys[temp] G += text(2,(xs[temp],f(x=xs[temp])+0.1),color='green') end S = (S*delx/3).n() show(G) exact = integrate(f,x,aa,bb) error = abs(exact-S).n() pretty_print(html("$\int f(x) dx =$"+str(exact.n()))) pretty_print("The Simpson's approximation to the integral using "+str(n)+" data points is "+str(S)+" and only has error "+ str(error))

$\int f(x) dx =$ 1.12176568093317

$\newcommand{\Bold}[1]{\mathbf{#1}}\verb|The|\phantom{\verb!x!}\verb|Simpson's|\phantom{\verb!x!}\verb|approximation|\phantom{\verb!x!}\verb|to|\phantom{\verb!x!}\verb|the|\phantom{\verb!x!}\verb|integral|\phantom{\verb!x!}\verb|using|\phantom{\verb!x!}\verb|24|\phantom{\verb!x!}\verb|data|\phantom{\verb!x!}\verb|points|\phantom{\verb!x!}\verb|is|\phantom{\verb!x!}\verb|1.12176865946635|\phantom{\verb!x!}\verb|and|\phantom{\verb!x!}\verb|only|\phantom{\verb!x!}\verb|has|\phantom{\verb!x!}\verb|error|\phantom{\verb!x!}\verb|2.97853317654528e-6|$

1.12176568093317

Let's consider Gaussian Quadrature

var('c1,c2,w1,w2') eq1 = c1+c2 == 2 eq2 = c1*w1+c2*w2 == 0 eq3 = c1*w1^2+c2*w2^2 == 2/3 eq4 = c1*w1^3+c2*w2^3 == 0 solve([eq1,eq2,eq3,eq4],(c1,c2,w1,w2))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left[c_{1} = 1, c_{2} = 1, w_{1} = -\frac{1}{3} \, \sqrt{3}, w_{2} = \frac{1}{3} \, \sqrt{3}\right], \left[c_{1} = 1, c_{2} = 1, w_{1} = \frac{1}{3} \, \sqrt{3}, w_{2} = -\frac{1}{3} \, \sqrt{3}\right]\right]$

f = 1/(1+x) a = 0 b = 4 w1 = -1/sqrt(3) w2 = 1/sqrt(3) x1 = (b-a)*(w1+1)/2 x2 = (b-a)*(w2+1)/2 exact = integral(f,(x,a,b)) pretty_print(n(exact)) GQ = (f(x=x1) + f(x=x2))*(b-a)/2 pretty_print(n(GQ))

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

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

var('c1,c2,c3,w1,w2,w3') eq1 = c1+c2+c3 == 2 eq2 = c1*w1+c2*w2+c3*w3 == 0 eq3 = c1*w1^2+c2*w2^2+c3*w3^2 == 2/3 eq4 = c1*w1^3+c2*w2^3+c3*w3^3 == 0 eq5 = c1*w1^4+c2*w2^4+c3*w3^4 == 2/5 eq6 = c1*w1^5+c2*w2^5+c3*w3^5 == 0 solve([eq1,eq2,eq3,eq4,eq5,eq6,w2],(c1,c2,c3,w1,w2,w3))

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left[c_{1} = \left(\frac{5}{9}\right), c_{2} = \left(\frac{8}{9}\right), c_{3} = \left(\frac{5}{9}\right), w_{1} = -\frac{1}{5} \, \sqrt{5} \sqrt{3}, w_{2} = 0, w_{3} = \frac{1}{5} \, \sqrt{5} \sqrt{3}\right], \left[c_{1} = \left(\frac{5}{9}\right), c_{2} = \left(\frac{8}{9}\right), c_{3} = \left(\frac{5}{9}\right), w_{1} = \frac{1}{5} \, \sqrt{5} \sqrt{3}, w_{2} = 0, w_{3} = -\frac{1}{5} \, \sqrt{5} \sqrt{3}\right]\right]$

f = 1/(1+x) a = 0 b = 4 w1 = -sqrt(3/5) w2 = 0 w3 = sqrt(3/5) x1 = (b-a)*(w1+1)/2 x2 = (b-a)*(w2+1)/2 x3 = (b-a)*(w3+1)/2 exact = integral(f,(x,a,b)) pretty_print(n(exact)) GQ = (5/9*f(x=x1) + 8/9*f(x=x2) + 5/9*f(x=x3))*(b-a)/2 pretty_print(n(GQ))

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

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

# Do an attempt at an indefinite integral .... using trapezoidal rule # Let's start over with function stuff var('x') f = 1/(1+x) a = 2 b = 20 C = 0 # This is a guess exact = ln(1+x)+C n = 10 delx = (b-a)/n xs = [a] I = [C] pts = [(a,C)] for k in range(n): xs.append(xs[k]+delx) show(((f(x = xs[k]) + f(x = xs[k+1])))*delx/2) trap = (f(x = xs[k]) + f(x = xs[k+1]))*delx/2 I.append(I[k]+trap) pts.append((xs[k+1],I[k+1])) end show(I) S = spline(pts) exact = ln(1+x)-ln(3) (plot(S,x,a,b)+plot(exact,x,a,b,color="green")).show()

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{39}{80}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{57}{176}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{75}{308}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{93}{476}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{111}{680}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{129}{920}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{147}{1196}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{165}{1508}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{183}{1856}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{201}{2240}$

$\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, \frac{39}{80}, \frac{357}{440}, \frac{3249}{3080}, \frac{65463}{52360}, \frac{7401}{5236}, \frac{1871091}{1204280}, \frac{26248413}{15655640}, \frac{810880527}{454013560}, \frac{6845166981}{3632108480}, \frac{3585541929}{1816054240}\right]$

integrate(1/(1+x),x,2,3).n()

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

(7/24).n()

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

381 - A11 - Numerical Integration

365 days ago by Professor381