Computer Project 3: Numerical Integration
Note: Make sure to use the matlab ’format long e’ to print you
esults. Write a report that contains for each problem the matla
program, results, and discussion of results
Problem I:
1. Write a matlab m-file program trapez.m as trapez(’f’,a,b,n) and returns
the trapezoidal approximation T (f, h), h = (b−a)/n. Define your func-
tion in a m-file as
function f=myf(x)
f=....
Use n = 15 and run trapez to compute the integral:∫ 1
−0.2
exp(x2)dx
Print the result.
2. Write a matlab m-file program romberg.m which should be executed by
calling romberg(’f’,a,b,delta,nmax) which in turn uses trapez.m, with
n = 2m, to compute Romberg approximations. Your program should
stop when either |R(m,m) − R(m − 1,m − 1)| < delta or m > nmax
and return the approximation R(m,m). Note: your function should be
defined in a separate m-file f.m
function f = f(x)
f=..........
Inside trapez.m you should call to obtain f(x) as
feval(f,x)
Use delta = 10−8 and nmax = 30 and test you program on the integral∫ 3
1
cos(2x2)exp(x)dx
1
3. Print the Romberg table.
Problem II:
1. Write a function as an m.file mygauss.m that should be called as my-
gauss(’f’,a,b,ng) to approximate the integrate
∫
a f(x)dx using ng Gauss
points. The mfile should start with:
function I=mygauss(f,a,b,ng)
...
...
I=
define ξ and their co
esponding weights w from the file gauss.m posted
in Scholar under Resources for ng=5. Shift them to [a,b] and use the
Gauss-Legendre quadrature. Hint: inside mygauss.m use feval(f,s) to
compute f(s) for a given value s.
2. Test your program on I1 =
∫ 1
0 x
4dx to obtain a zero e
or.
3. Run your program to compute the following integral and print you
esult
I2 =
∫ 2
1
5cos(2x2)exp(x)dx.
4. Consider the curve described by the parametric equations (x(t) =
exp(t)cos(t2), y(t) = sin(5t)), t ∈ [0, 2] and define its length as l =∫ 2
0
√
x′2 + y′2dt. Divide the interval [0, 2] into m = 2, 4, 8, 16, 32, 64 uni-
form subintervals and use the composite Gauss quadrature with ng=4
to compute the length of the curve. Compute the e
ors for each case
y assuming the true length using m = 50, ng = 10. Print a table
with the first column contains the values of m, the second contains the
lengths of the curve and the third column contains the absolute value
of the e
ors.
5. Print and discuss all results.
Problem III:
2
1. Write a function as an m.file mygauss2d.m that should be called as my-
gauss2d(’f2d’,ng) to approximate the double integral
∫
a
∫ d
c f(x, y)dydx.
on the domain 0 < x < 2.5, cos(x) < y < 2ex using ng × ng Gauss
points.
The mfile mygauss2d.m should start with:
function I=mygauss2D(f,ng)
...
I=
2. define an mfile f2d.m such that
function f = f2d(x,y)
...
to define the function f(x, y) = (x2 + y2)cos(x)
3. run your program and print the double integral of f using 2 by 2 and
5 by 5 and 10 by 10 Gauss points. Discuss your results.
4. Composite Gauss: (Bonus Question=+5)
Divide the domain of integration into four regions using the vertical line
x = (a + b)/2 and the curve y = (c(x) + d(x))/2 apply composite 2d
Gauss integration with 10 points to compute the double integral given
above. You may modify your function to mygauss2d(’f2d’,a,b,’c’,’d’,ng)
where c and d are defined a functions in m-files c.m and d.m and ng is
the number of Gauss points in each variable.
Note: Your are not allowed to discuss programming issues with
anyone except your instructor.
3