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 GaussLegendre 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 GaussLaguerre form put
t = 3(x  2), dt = 3dx
to get
(c) To put the integral in GaussChebyshev form put cos u = x, sin u du = dx and use the results that cos(2u) = 2cos^{2 }u  1 = 2x^{2}  1 and
^{}
3. (a) Since tan x = x + O(x^{3}) the leading term in the integrand is x^{1p} and the integral exists for
p < 2. (Ignore the singularity at x = p/2).
(b) Since 1  cos x = x^{2} + O(x^{4}) the leading term in the integrand is x^{2p} and the integral exists for p < 3.
(c) Here the singularity is at x = 2. Factoring the numerator the integrand becomes
(2  x)^{1/2p}(2 + x)^{1/2 }and the integral exists for p < 3/2.
4.(a) The integrand behaves as x^{2q} for large x so that the integral exists for 2  q < 1 or q > 3.
(b) The integrand behaves as e^{(4q)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 e^{4x}.
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(x^{16/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 gausslegendre integration
c
implicit real*8(ah,oz)
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(ah,oz)
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(ah,oz)
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 gausslaguerre integration
c
implicit real*8(ah,oz)
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(ah,oz)
e = dexp(1.d0)
func = (x+6.0d0)**3/8.1d1/e**6
return
end
integral is 0.11200288D01
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(ah,oz)
external func
open(unit=2,file='genint.out',status='new')
errabs = 0.d0
errrel = 1.d6
irule = 2
write(*,1)
1 format( 'input limits of integration')
read(*,*)a,b
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(ah,oz)
func = dsin(x)/(1.d0+dcos(x))
return
end
integral has the value 0.548230D+00
estimated error: 0.609D14
8 (b) Use of IMSL routine DQDAGI
c general gaussian integration over infinite range
implicit real*8(ah,oz)
external func
open(unit=2,file='geninti.out',status='new')
errabs = 0.d0
errrel = 1.d6
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(ah,oz)
func = dexp(x*x)/(1.d0 + x*x)
return
end
integral from infinity to infinity has the value 0.134329D+01
estimated error: 0.741D06
Assignment 3  Problem 6  FW03
#include <math.h>
#include <imsl.h>
#define QUADPTS 10
double func (double x){
return (sqrt(x*x + 6.0*x + 11.0)/sqrt(8.0));
}
main()
{
int j;
double weights[QUADPTS], points[QUADPTS], sum;
/* Produce the Gauss Legendre quadrature points */
imsl_d_gauss_quad_rule (QUADPTS, weights, points, 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
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>
#define QUADPTS 10
double func (double x){
double e = exp(1.0);
return (pow((x+6.0),3.0)/81.0/pow(e,6));
}
main()
{
int j;
double weights[QUADPTS], points[QUADPTS], alpha, sum;
/* Produce the Gauss Legendre */
/* quadrature points */
alpha =0.0;
imsl_d_gauss_quad_rule (QUADPTS, weights, points,
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.048e313