% Script 2 for Week 8
t=linspace(.01,2);
g= @(x) exp(-x.^2).*log(x).^2
y=g(t);
plot(t,y)
integral(@(x) exp(-x.^2).*log(x).^2,0,Inf)
% Place the trapuneq.m M-file in your working directory:
I = trapuneq(t,y)
%**********************************************************
% Integrate 1./(sqrt(x+y).*(1+x+y).^2 over the triangle
% 0 <= x <= 1,
% 0 <= y <= 1-x.
% The integrand is infinite at (0,0).
% The true value of the integral is pi/4 - 1/2.
fun = @(x,y) 1./( sqrt(x + y) .* (1 + x + y).^2 );
ymax = @(x) 1 - x;
Q = integral2(fun,0,1,0,ymax)
pi/4 - 1/2
%**********************************************************
% Integrate over the unit sphere
f3 = @(x,y,z) x.*cos(y) + x.^2.*cos(z)
xmin = -1;
xmax = 1;
ymin = @(x)-sqrt(1 - x.^2);
ymax = @(x) sqrt(1 - x.^2);
zmin = @(x,y)-sqrt(1 - x.^2 - y.^2);
zmax = @(x,y) sqrt(1 - x.^2 - y.^2);
Q = integral3(f3,xmin,xmax,ymin,ymax,zmin,zmax)
%**********************************************************
% trapuneq works the same way as the built-in MATLAB function trapz but,
% of course, trapz has more options. I give trapuneq so you can see how
% trapz works in default configuration.
% trapz(t,y) gives the same answer as trapuneq(t,y).
% trapz(y) with y by itself presumes that the domain increments by 1 from
% y value to y value. If the ACTUAL increment is constant increment d then
% the approximate integral would be d*trapz(y) which is the same as trapz(d*y).
% Given independent data
t=[0 0.1 0.3 0.5 0.7 0.8 0.9 1.1 1.3 1.4 1.5 1.6 1.8 1.9]
% integrate the dependent data
f=@(x) sin(x)+cos(x)
z=f(t);
% over the interval.
tun=trapuneq(t,z)
exact=integral(@(x) sin(x)+cos(x),0,1.9)
x=linspace(0,1.9,114);
y=f(x);
tr=trapz(x,y)
% The 1/3 Simpson rule is 5th order accurate but requires an odd number of data
% points. If you have an even number of data points you can use trapz on
% the first segment or ThreeEighthSimp on the first 3 segments and patch
% the integral together (add the pieces) using OneThirdSimp on the
% (now even number of) remaining data points. Note the effect that just
% one segment of 3rd order accurate trapz has on the error calculated below.
% 3/8 Simpson is 5th order accurate and should always be used instead
% for that initial segment.
% Though both the provided 3/8 Simpson and 1/3 Simpson functions have been
% designed to work with any number of data points, we illustrate the
% importance of paying attention to effects like this below.
ts=trapz(x(1:2),y(1:2))+simp13(x(2:114),y(2:114))
sim=simp38(x(1:4),y(1:4))+simp13(x(4:114),y(4:114))
tun2=trapuneq(x,y)
(tun2-exact)/exact
(tr-exact)/exact
(ts-exact)/exact
(sim-exact)/exact
%**********************************************************
% There also is the equivalent of an antiderivative, given by cumtrapz.
% cumtrapz(t,y) is a trapezoid-rule approximation to a function whose
% derivative at the t values is the y values.
f=@(x) sin(x)+cos(x)
t=[0 0.1 0.3 0.5 0.7 0.8 0.9 1.1 1.3 1.4 1.5 1.6 1.8 1.9];
z=f(t);
% Following is the actual antiderivative of the function f that passes through
% the origin. cumtrapz(t,z) will always assign zero to the start of
% the t vector so may differ from another antiderivative by an
% additive constant..
antideriv=@(x)sin(x)-cos(x)+1
yyy=antideriv(t);
yy=cumtrapz(t,z);
plot(t,yyy,t,yy)
%**********************************************************
% Finally, I would like to expand on the MATLAB interpolation functions
% interp1, interp2 and interp3 for 1, 2 and 3 dimensional interpolation.
% Vq = interp1(X,V,Xq,METHOD) specifies the interpolation method.
% The available methods are:
% 'linear' - (default) linear interpolation
% 'nearest' - nearest neighbor interpolation
% 'next' - next neighbor interpolation
% 'previous' - previous neighbor interpolation
% 'spline' - piecewise cubic spline interpolation (SPLINE)
% 'pchip' - shape-preserving piecewise cubic interpolation
% 'cubic' - same as 'pchip'
% 'v5cubic' - the cubic interpolation from MATLAB 5, which does not
% extrapolate and uses 'spline' if X is not equally spaced.
% 'makima' - modified Akima cubic interpolation
% Vq = interp2(XX, YY, V, Xq, Yq, 'method')
% Comparing pchip with SPLINE:
% The function s(x) supplied by SPLINE is constructed in exactly the same way,
% except that the slopes at the X(j) are chosen differently, namely to make
% even D^2s(x) continuous. This has the following effects.
% SPLINE is smoother, i.e., D^2s(x) is continuous.
% SPLINE is more accurate if the data are values of a smooth function.
% pchip has no overshoots and less oscillation if the data are not smooth.
% pchip is less expensive to set up.
% which returns interpolated values of a function of two variables at
% specific query points using linear interpolation as default (no method
% specified). The results always pass through the original sampling of
% the function. XX and YY contain grids of X coordinates and Y coordinates,
% respectively, of the sample points. V contains the corresponding function
% values at each sample point. Xq and Yq contain the coordinates of the
% query points. They can be individual values or vectors of values to be
% interpolated.
% 'method' for BOTH interp2 and interp3 is
% 'nearest' - nearest neighbor interpolation
% 'linear' - bilinear interpolation
% 'spline' - spline interpolation
% 'cubic' - bicubic interpolation as long as the data is uniformly spaced
% 'makima' - modified Akima cubic interpolation
% Let's look at an example.
X=[0 5 10 15 20 25 30]
Y=[0 10 20]
[XX YY]=meshgrid(X,Y)
Table=[ 14.6000 12.8000 11.3000 10.1000 9.0900 8.2600 7.5600;
12.9000 11.3000 10.1000 9.0300 8.1700 7.4600 6.8500;
11.4000 10.3000 8.9600 8.0800 7.3500 6.7300 6.2000]
surf(XX,YY,Table)
V12_10_cubic = interp2(XX, YY, Table, 12, 10,'cubic')
% V12_10_cubic = 9.6526
Vlin = interp2(XX, YY, Table, 12, 15)
% Vlin = 9.1400
interp2(XX, YY, Table, 12, 15,'linear')
% ans = 9.1400
Vspl=interp2(XX, YY, Table, 12, 15,'spline')
% Vspl = 9.1024
Vcubic=interp2(XX, YY, Table, 12, 15,'cubic')
% Vcubic = 9.1036