>

restart:

 

M242  Lab4      

>	Digits:=18:  # computations done with more digits.

 

We'll do several things here. They are on the numerical methods. 

A. Error Analysis.  For Simpson's , the error is bounded by (K*(b-a)*h^4)/180, where K is the max of the fourth derivative of f(x) on [a,b]. 

As an example, let's look at the following integral (which is similar to #22 in 7.7): 

>	Int(exp(-x^2), x=0..1); value(%); bell:=evalf(%);

 

 

 

	(1)

 

The graph is the famous bell-shaped curve plotted on [-4, 4]. 

>	plot( exp(-x^2), x=-4..4);

 

If we use Simpson's with n = 4,  a lazy approach would consist of the following steps: 

>	with(Student[Calculus1]):

 

>	s4:= ApproximateInt(exp(-x^2), x=0..1, method = simpson, partition= 4):

 

>	evalf([s4, bell, abs(s4-bell)]);

 

(2)

 

If we did not know the exact value of bell, the error in using s4 to do the approximation would be estimated by looking at the maximum of the fourth derivative of the integrand on the given interval. 

We recall the error for Simpson's is given by  f""(c)*(b-a)^5/(180*n^4) where c is some number in (a,b). 

>	BS:= exp(-x^2):   # BS is an expression, not a function.

 

>	diff(BS,x$4); D4BS:= simplify(%);   # D4BS denotes the 4th derivative of B.

 

 

	(3)

 

>	plot( D4BS, x=0..1);

 

We see the maximum is 12 occurring at x = 0. Therefore the error estimate is 12/(180*4^4).  

>	e4:= evalf(12/180/4^4);

 

(4)

 

>	evalf([s4-e4, s4+e4]); bell; # the interval is the one where the exact value would be in.

 

 

	(5)

 

PROBLEM ONE: 

Repeat the procedures above  with n = 8 for another 'bell-shaped' function  4/(1 + x^2). 

. 

>	Int( 4/(1+x^2), x=0..1); exactvalue:= value(%);

 

 

	(6)

 

>

 

INTRODUCTION 

On a Newton-Cotes Method of integration. 

For Simpson's rule, what we did was to approximate the graph of a function over two consective intervals by a parabola and came up with the quadrature rule with weights [1, 4, 1]. 

Area := h/3*(y0 + 4*y1 + y2) .    See p.502. We also did it in class. 

Let's verify the quadrature rule with weights being [1, 3 , 3, 1] if we approximate the given function by a cubic polynomial  p:= A*x^3 + B*x^2 + C*x + D.  Since there are four unknowns, we need four points, i.e. three consective intervals. 

The new formula we want to verify is that Area:= 3*h/8*(y0 + 3*y1 + 3*y2 + y3) 

>

 

>	p:=x-> A*x^3 + B*x^2 + C*x + D;

 

(7)

 

>	int(p(x), x = 0..3*h);

 

(8)

 

>	3*h/8*(p(0) + 3*p(h) + 3*p(2*h) + p(3*h)); expand(%);

 

 

	(9)

 

So this quadrature formula suggests the following rule. Break up an interval [a, b] into 3*N number 

of subintervals. Evaluate a given function  h(x) at  a=x0, x1, x2, ..., b=x_3N . 

Form the sum with h:= (b - a)/(3*N) . 

3*h/8*( y0 + 3*y1 + 3*y2 + 2*y3 + 3*y4 + 3*y5 + 2*y6 + 3*y7 + ... + 3*y_(3N-1) + y_(3N) ). 

Let's apply this formula to the famous elliptic integral: 

>	h:=x-> 1/sqrt(9*cos(x)^2 + 16*sin(x)^2); ell:= int(h(x), x=0..Pi/2);

 

 

	(10)

 

PROBLEM TWO: 

Write the sum above as a procedure in terms of N, a and b.  Evaluate the sum for (i) N = 4 (which corresponds to 12 sub-intervals) and (ii) N = 5 , i.e.15 intervals. Compare the exact value which is 'ell' with your sum each time. 

For your reference, we saw this proc in the last lab. We also adapted the proc to the Simpson's rulr for practice. 

----------------------------------------------------------------------------- 

>	F:=x-> 4/(1+x^2);

 

(11)

 

>	T:= proc(N, a, b) local h, k, s, exva; h:= (b-a)/N: s:= 2*sum(f(a+k*h), k=1..N-1): exva:= int( F(x), x=a..b): print( exva,  evalf([exva, (F(a) + s + F(b))*h/2])); end:

 

>	T(100,0,1);

 

(12)

 

>

 

>