MATH3242.03/COSC3122.03

FW99

Solutions to Assignment 3

1. To evaluate the integral correct to 4 decimal places (i.e. with an error e < 5´10-5) we use adoptive integration with Simpson's rule. If we denote the elementary Simpson's rule over the interval [a, b] by S(a, b) we get the following table of values. The error is calculated from the formula [S(a,b) - S(a,(a+b)/2) - S((a+b)/2,b)]/15. The entries denoted by * have reached the required accuracy.

 Error S(1, 2) 0.532827 S(1, 3/2) 0.468284 0.496393 { 0.002430 S(3/2, 2) 0.028109 S(1, 5/4) 0.235692 0.467992* { 0.000020 S(5/4, 3/2) 0.232300 S(3/2, 7/4) 0.115758 0.026621 { 0.000100 S(7/4, 2) -0.089136 S(3/2, 13/8) 0.079912 0.115707* { <0.00001 S(13/8, 7/4) 0.035795 S(7/4, 15/8) -0.017888 -0.089163* { <0.00001 S(15/8, 2) -0.071275

Result for integral (sum of * values) = 0.494536.

2(a). To put the first integral in standard Gauss-Legendre form make the substitution

x = (t + 3)/2, dx = dt/2

to get

For the second integral put

t = 2u+2, dt = 2du

to get

(b) To put the integral in Gauss-Laguerre form put

t = 3(x - 2), dt = 3dx

to get

(c) To put the integral in Gauss-Chebyshev form put cos u = x, sin u du = dx and use the results that cos(2u) = 2cos2 u - 1 = 2x2 - 1 and

3. (a) Since tan x = x + O(x3) the leading term in the integrand is x1-p and the integral exists for

p < 2. (Ignore the singularity at x = p/2).

(b) Since 1 - cos x = x2 + O(x4) the leading term in the integrand is x2-p and the integral exists for p < 3.

(c) Here the singularity is at x = 2. Factoring the numerator the integrand becomes

(2 - x)1/2-p(2 + x)1/2 and the integral exists for p < 3/2.

4.(a) The integrand behaves as x2-q for large x so that the integral exists for 2 - q < -1 or q > 3.

(b) The integrand behaves as e(4-q)x for large x if q > 0 so the integral exists for 4 - q < 0 or q > 4.

The integral will not exist for q < 0 since then the integrand behaves as e4x.

5 Since the error in Simpson's rule involves the fourth derivative we must treat the singular part using a polynomial of degree n > 4 + 14/3 -1 = 8

The Taylor polynomial for (1 - cos x)2 is

The second integral is evaluated by the elementary Simpson's Rule yielding -0.000023. Note that the integrand is O(x16/3) and hence is zero at x = 0. The original integral then has the approximate value 0.732841.

c Assignment 3 - question 6 fw03

c standard gauss-legendre integration

c

implicit real*8(a-h,o-z)

dimension qx(10), qw(10), qfix(1)

open(unit=2,file='gleg.out',status='new')

n = 5

alpha = 0.d0

beta = 0.d0

nfix = 0

iweigh = 1

call dgqrul(n, iweigh, alpha, beta, nfix, qfix, qx, qw)

sum = 0.d0

do j = 1,n

sum = sum + qw(j)*func(qx(j))

end do

write(2,3) sum, n

3 format(' integral is',d16.8,' with',i3,' points')

stop

end

c

c integrand for first integral in question 2(a)

function func(x)

implicit real*8(a-h,o-z)

func = dsqrt(x*x+6.d0*x+1.1d1)/dsqrt(8.d0)

return

end

integral is 0.23519513D+01 with 5 points

integral is 0.23519513D+01 with 10 points

c integrand for second integral in question 2(a)

function func(x)

implicit real*8(a-h,o-z)

func = 2.d0*dtan(6.d0*x+5.d0)

return

end

integral is -0.10489903D+02 with 5 points

integral is 0.35649177D+01 with 10 points

These results do not seem to be consistent. This should alert you to the fact that there is something wrong. In fact, the integrand has several singularities in the range of integration and thus the integral is not defined.

c Assignment 3 - question 7- fw03

c standard gauss-laguerre integration

c

implicit real*8(a-h,o-z)

dimension qx(10), qw(10), qfix(1)

open(unit=2,file='glag.out',status='new')

n = 10

alpha = 0.d0

beta = 0.d0

nfix = 0

iweigh = 6

call dgqrul(n, iweigh, alpha, beta, nfix, qfix, qx, qw)

sum = 0.d0

do j = 1,n

sum = sum + qw(j)*func(qx(j))

end do

write(2,3) sum

3 format(' integral is',d16.8)

stop

end

c

c integrand for question 2b)

function func(x)

implicit real*8(a-h,o-z)

e = dexp(1.d0)

func = (x+6.0d0)**3/8.1d1/e**6

return

end

integral is 0.11200288D-01

N.B. This is the exact result since the integrand excluding the weight function is a polynomial of degree 3.

8 (a) Use of IMSL routine DQDAG

c general gaussian integration over finite range

implicit real*8(a-h,o-z)

external func

open(unit=2,file='genint.out',status='new')

errabs = 0.d0

errrel = 1.d-6

irule = 2

write(*,1)

1 format( 'input limits of integration')

call dqdag(func,a,b,errabs,errrel,irule,result,errest)

write(2,2) result, errest

2 format(' integral has the value',d16.6/' estimated error:',d10.3)

stop

end

c integrand - ass3 - q8(a)

function func(x)

implicit real*8(a-h,o-z)

func = dsin(x)/(1.d0+dcos(x))

return

end

integral has the value 0.548230D+00

estimated error: 0.609D-14

8 (b) Use of IMSL routine DQDAGI

c general gaussian integration over infinite range

implicit real*8(a-h,o-z)

external func

open(unit=2,file='geninti.out',status='new')

errabs = 0.d0

errrel = 1.d-6

interv = 2

bound = 0.

call dqdagi(func,bound,interv,errabs,errrel,result,errest)

write(2,2) result, errest

2 format(' integral from -infinity to infinity has the value',

1 d16.6/' estimated error:',d10.3)

stop

end

c integrand - ass3 - q9b

function func(x)

implicit real*8(a-h,o-z)

func = dexp(-x*x)/(1.d0 + x*x)

return

end

integral from -infinity to infinity has the value 0.134329D+01

estimated error: 0.741D-06

Assignment 3 - Problem 6 - FW03

#include <math.h>

#include <imsl.h>

double func (double x){

return (sqrt(x*x + 6.0*x + 11.0)/sqrt(8.0));

}

main()

{

int j;

/* Produce the Gauss Legendre quadrature points */

/* integrate the function func */

sum = 0.0;

for(j = 0; j < QUADPTS; j++) {

sum += weights[j]*func(points[j]);

}

printf("The integral from -1 to 1 of func is\n");

printf(" %10.6f with %d points\n", sum, QUADPTS);

}

The integral from -1 to 1 of func is

2.351951 with 5 points

The integral from -1 to 1 of func is

2.351951 with 10 points

Function for second integral

double func ( double x){

return(2.0*tan(6.0*x+5.0));

}

The integral from -1 to 1 of func is

-10.489903 with 5 points

The integral from -1 to 1 of func is

3.564918 with 5 points

Note inconsistent results due to singularities in integrand.

Assignment 3 - Problem 7 - FW03

#include <math.h>

#include <imsl.h>

double func (double x){

double e = exp(1.0);

return (pow((x+6.0),3.0)/81.0/pow(e,6));

}

main()

{

int j;

/* Produce the Gauss Legendre */

alpha =0.0;

IMSL_GEN_LAGUERRE, alpha, 0);

/* integrate the function func */

sum = 0.0;

for(j = 0; j < QUADPTS; j++) {

sum += weights[j]*func(points[j]);

}

printf("The integral from -1 to 1 of func is\n");

printf(" %10.6f with %d points\n", sum, QUADPTS);

}

The integral from -1 to 1 of func is

0.011200 with 10 points

Assignment 3 - Problem 8a - FW03

#include <math.h>

#include <imsl.h>

double fcn(double x);

double q, errest;

main()

{

/* evaluate the integral */

q = imsl_d_int_fcn (fcn, 0.0, sqrt(2.0), IMSL_ERR_EST, errest, 0);

/* print the result */

printf("integral = %10.6f with estimated error %10.3e\n", q, errest);

}

double fcn(double x)

{

float y;

y = sin(x)/(1.0 + cos(x));

return y;

}

integral = 0.548230 with estimated error 0.000e+00

Assignment 3 - Problem 8b _ FW03

#include <math.h>

#include <imsl.h>

double fcn(double x);

main()

{

double q, errest;

/* Evaluate the integral */

q = imsl_d_int_fcn_inf (fcn, 0.0,

IMSL_INF_INF,0.0,

IMSL_ERR_EST, errest, 0);

/* Print the result */

printf("integral = %10.6f with estimated error %10.3e\n", q, errest);

}

double fcn(double x)

{

return (exp(-x*x)/(1.0 + x*x));

}

integral = 1.343293 with estimated error 1.048e-313