% 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