% started off trying out matlab's numerical quadrature commands quad and quad8 % Both use adaptive quadrature >> help quad QUAD Numerically evaluate integral, low order method. Q = QUAD('F',A,B) approximates the integral of F(X) from A to B to within a relative error of 1e-3 using an adaptive recursive Simpson's rule. 'F' is a string containing the name of the function. Function F must return a vector of output values if given a vector of input values. Q = Inf is returned if an excessive recursion level is reached, indicating a possibly singular integral. Q = QUAD('F',A,B,TOL) integrates to a relative error of TOL. Use a two element tolerance, TOL = [rel_tol abs_tol], to specify a combination of relative and absolute error. Q = QUAD('F',A,B,TOL,TRACE) integrates to a relative error of TOL and for non-zero TRACE traces the function evaluations with a point plot of the integrand. Q = QUAD('F',A,B,TOL,TRACE,P1,P2,...) allows parameters P1, P2, ... to be passed directly to function F: G = F(X,P1,P2,...). To use default values for TOL or TRACE, you may pass in the empty matrix ([]). See also QUAD8, DBLQUAD. % note the possibility of getting a trace, which, in this case, means, % a picture showing the actual points on the function used by the command >> f = inline('sqrt(x)') f = Inline function: f(x) = sqrt(x) % integrate the function over the interval [0..1]. Note that the integrand % has infinite slope at the left endpoint, making it difficult to integrate % numerically. >> quad(f,0,1,1e-6,1) Warning: Recursion level limit reached in quad. Singularity likely. > In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 102 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m at line 72 Warning: Recursion level limit reached 40 times. > In C:\MATLABR11\toolbox\matlab\funfun\quad.m at line 80 ans = 0.6667 % let's let up on the tolerance a bit, using 1e-3 rather than 1e-6: >> quad(f,0,1,1e-3,1) Warning: Recursion level limit reached in quad. Singularity likely. > In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 102 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m (quadstp) at line 141 In C:\MATLABR11\toolbox\matlab\funfun\quad.m at line 72 Warning: Recursion level limit reached 8 times. > In C:\MATLABR11\toolbox\matlab\funfun\quad.m at line 80 ans = 0.6667 % still difficult. Look at the answer in more detail: >> format long ans ans = 0.66666250953207 % the graph shows high density of points near the left endpoint. % Now try out quad8 which is supposedly more efficient. >> quad8(f,0,1,1e-3,1) Warning: Recursion level limit reached in quad8. Singularity likely. > In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 129 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m at line 70 Warning: Recursion level limit reached 2 times. > In C:\MATLABR11\toolbox\matlab\funfun\quad8.m at line 77 ans = 0.66666663842017 % still difficulty, but the answer is slightly better; in fact, way better % than the tolerance we specified. This may point to a not very sophisticated % assignment of the tolerance when subintervals are considered. % notice the points are not uniformly spaced, in fact, there % are strange gaps. % These are due to use of Gauss rules which are open rules. In fact, from % the graph, we guess that 3-point Gauss-Legendre rules are being used. % Let's try an even less smooth function, namely a step function: >> f = inline('x>.3') f = Inline function: f(x) = x>.3 >> quad8(f,0,1,1e-3,1) Warning: Recursion level limit reached in quad8. Singularity likely. > In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 129 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 124 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 124 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 124 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 124 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m (quad8stp) at line 123 In C:\MATLABR11\toolbox\matlab\funfun\quad8.m at line 70 Warning: Recursion level limit reached 2 times. > In C:\MATLABR11\toolbox\matlab\funfun\quad8.m at line 77 ans = 0.69991207492835 % the integral should be .7, so, this is pretty good! % The book also considers integration of spline functions, but does it wrong, % I think. Better to do the following: Suppose we have in hand the command % fnint(pp) that returns to us the piecewise polynomial function that is the % integral or primitive of the piecewise polynomial function f in pp. % Then we can compute integral f(x) dx over [a..b] in one statement: % diff( ppval( fnint(pp) , [a b] )) % Let's try that out: >> x = linspace(0,2*pi,5); >> sp = spline(x,sin(x)); >> diff(ppval(fnint(sp),[0,2*pi])) ans = 0 % That's remarkable! The exact answer! (Actually, it's just symmetry.) % Let's try it instead on the interval >> diff(ppval(fnint(sp),[0,pi])) ans = 2.09439510239320 % .. while the exact value should be 2 . % % Check this against the spline interpolant to 5 points on the interval % [0 .. pi]. That makes the breakpoint space half of what it was, hence, % since the error in spline interpolation behaves like spacing^4, this should % cut the error by 2^4 = 16. >> x = linspace(0,pi,5); >> sp = spline(x,sin(x)); >> diff(ppval(fnint(sp),[0,pi])) ans = 2.00455975498442 % compare this with 1/16 of the earlier error: >> .094395/16 ans = 0.00589968750000 % ... pretty close! %% the rest of the class time was spent on matrices in matlab, especially % matrix multiplication. See the notes under the story so far. % get some matrices to play with: >> a = reshape(1:10,2,5) a = 1 3 5 7 9 2 4 6 8 10 >> b = reshape(1:6,3,2) b = 1 4 2 5 3 6 % well, to multiply a with b, I need a to have only three columns. % Here is one way to get rid of the last two columns: >> a(:,[4 5]) = [] a = 1 3 5 2 4 6 % now I compute the outer product of the first column of a with the first % row of b , to make clear that I don't get a scalar but a full matrix: >> a(:,1)*b(1,:) ans = 1 4 2 8 % Of course, if I multiply the row with the column, i.e., the other way around % I get the scalar or inner product: >> b(1,:)*a(:,1) ans = 9 %% Then we generated the Hilbert matrix, the one with A(i,j) = 1/(i+j-1). %% We don't want to use any loops !!!! function H = hilb(n) %HILB Hilbert matrix. % HILB(N) is the N by N matrix with elements 1/(i+j-1), % which is a famous example of a badly conditioned matrix. % See INVHILB for the exact inverse. % % This is also a good example of efficient MATLAB programming % style where conventional FOR or DO loops are replaced by % vectorized statements. This approach is faster, but uses % more storage. % C. Moler, 6-22-91. % Copyright (c) 1984-98 by The MathWorks, Inc. % $Revision: 5.6 $ $Date: 1997/11/21 23:28:56 $ % I, J and E are matrices whose (i,j)-th element % is i, j and 1 respectively. J = 1:n; J = J(ones(n,1),:); I = J'; E = ones(n,n); H = E./(I+J-1); %%%%% this can be done slightly more swiftly now. E.g., for n=3: >> n = 3; >> J = repmat(1:n,n,1); >> H = 1./(J'+J-1) H = 1.00000000000000 0.50000000000000 0.33333333333333 0.50000000000000 0.33333333333333 0.25000000000000 0.33333333333333 0.25000000000000 0.20000000000000 % I contacted Cleve Moler about this. His present solution is: >> [I,J] = ndgrid(1:n); >> H = 1./(I+J-1) % ... very neat! diary off