381 - A11 - Numerical Integration

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

                                
                            

                                

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()) 
       
1.70710678118655
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.&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)) 
       
1.12176568093317

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

                                
                            

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

                                
                            

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

                                
                            

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

                                
                            

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

                                
                            

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

                                
                            

                                
(7/24).n()