# 381 - A11 - Numerical Integration

## 282 days ago by Professor381

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

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

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)

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 1.707106781186551.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)

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.707106781186551.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.&nbsp; Breaking up a larger interval therefore into pieces and reusing each formula over and over on the smaller pieces might be advantageous.&nbsp; 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 1.12176568093317