< M A T L A B > Copyright 1984-2001 The MathWorks, Inc. Version 6.1.0.450 Release 12.1 May 18 2001 To get started, type one of these: helpwin, helpdesk, or demo. For product information, visit www.mathworks.com. % Following indicates what happens if you don't enter the syms command >> f = inline(vectorize((1-cos(y-y^2))/y),'y') ??? Undefined function or variable 'y'. >> syms u v x y z >> f = inline(vectorize(sin(x^3 - 7*x)),'x') f = Inline function: f(x) = sin(x.^3-7.*x) >> quad(f,2,4) ans = 0.1734 >> quadl(f,2,4) ans = 0.1734 % The answers above agree as they should. % We can be pretty certain the integral is about .1734 >> f = inline(vectorize(sqrt(x^2 + 4/y^2)),'x','y') f = Inline function: f(x,y) = (x.^2+4./y.^2).^(1./2) >> dblquad(f,3,4,1,2) ans = 3.7749 % Now let's use method 1, a change of variables to a rectangular region >> g = inline(vectorize((1-u)*sin(u*v*(1-u))),'u','v') g = Inline function: g(u,v) = (1-u).*sin(u.*v.*(1-u)) >> dblquad(g,0.1.0.1) ??? dblquad(g,0.1.0.1) | Error: ")" expected, "numeric value" found. % Ooops, I typed periods instead of commas. >> dblquad(g,0,1,0,1) ans = 0.0415 % Okay, that's better. Now let's do the inner integration by hand >> f = inline(vectorize((1-cos(y-y^2))/y),'y') f = Inline function: f(y) = (1-cos(-y+y.^2))./y >> quad(f,0,1) Warning: Divide by zero. > In /usr/local/matlab6.1/toolbox/matlab/funfun/inlineeval.m at line 13 In /usr/local/matlab6.1/toolbox/matlab/funfun/@inline/feval.m at line 34 In /usr/local/matlab6.1/toolbox/matlab/funfun/quad.m at line 62 ans = 0.0415 % So we get the same answer as before. The warnings tell you to be suspicious of your answer. % But we can be pretty confident of its accuracy since it agrees with the method 1 answer. % Now let's use ezint241 >> ezint241(sin(x*y), x,0,1-y, y,0,1) Integrate sin(x.*y) dxdy with x limits 0 to 1-y and y limits 0 to 1 ans = 0.0415 % and again we get the same answer. % Using method 2, doing theinner integration by hand we get: >> ezint241(cos(x^2*exp(y))-cos(x^2*exp(y)+6-3*x-2*y), y,0, 3-1.5*x, x,0,2) Integrate cos(x.^2.*exp(y))-cos(x.^2.*exp(y)+6-3.*x-2.*y) dydx with y limits 0 to 3-3./2.*x and x limits 0 to 2 ans = 2.4396 % Now do the triple integral directly >> ezint241(sin(x^2*exp(y)+z), z,0,6-3*x-2*y, y,0,3-1.5*x, x,0,2) Integrate sin(x.^2.*exp(y)+z) dzdydx with z limits 0 to 6-3.*x-2.*y with y limits 0 to 3-3./2.*x and x limits 0 to 2 ans = 2.4406 % The answers are close, the default tolerance is only .5% error for triple integrals % and our two answers are well within this tolerance % So we figure the integral is about 2.44 --------------------------- % Now for some extra stuff. % Here is an example of finding the area of a parameterized surface. % First tell matlab the parameterization, % In this example the surface is parameterized by % F(u,v) = uv i + e^{u+v} j + sin(u/v) k, 1<=u<=2, 3<=v<=4. >> F=[u*v, exp(u+v), sin(u/v)] F = [ u*v, exp(u+v), sin(u/v)] % Now calculate the partial derivatives of F >> Fu = diff(F,u) Fu = [ v, exp(u+v), cos(u/v)/v] >> Fv= diff(F,v) Fv = [ u, exp(u+v), -cos(u/v)*u/v^2] % Matlab can calculate the cross product >> cross(Fu,Fv) ans = [ -exp(u+v)*cos(u/v)*u/v^2-cos(u/v)/v*exp(u+v), 2*cos(u/v)/v*u, v*exp(u+v)-exp(u+v)*u] % I forgot to give a name to the cross product so I can refre to it later % The following fixes that, the variable ans is always the last answer >> K = ans K = [ -exp(u+v)*cos(u/v)*u/v^2-cos(u/v)/v*exp(u+v), 2*cos(u/v)/v*u, v*exp(u+v)-exp(u+v)*u] % Now I need to find the length, the square root of the dot product of K with itself % Matlab's function dot is meant to work for complex vectors also which makes it % look pretty complicated, with all those conj terms (that means take the complex conjugate, % which does nothing to real vectors. >> dot(K,K) ans = conj(-exp(u+v)*cos(u/v)*u/v^2-cos(u/v)/v*exp(u+v))*(-exp(u+v)*cos(u/v)*u/v^2-cos(u/v)/v*exp(u+v))+4*conj(cos(u/v)/v*u)*cos(u/v)/v*u+conj(v*exp(u+v)-exp(u+v)*u)*(v*exp(u+v)-exp(u+v)*u) % The following way of taking the dot product gives better looking answers % for real vectors, since it doesn't have all those conj terms % If you want to know what it does, K.*K squares all the entries of K, % and sum adds all the entries up. >> sum(K.*K) ans = (-exp(u+v)*cos(u/v)*u/v^2-cos(u/v)/v*exp(u+v))^2+4*cos(u/v)^2/v^2*u^2+(v*exp(u+v)-exp(u+v)*u)^2 % Now we can find the area which we get by integrating the length of K. >> ezint241(sqrt(sum(K.*K)), u, 1,2, v, 3,4) Integrate ((-exp(u+v).*cos(u./v).*u./v.^2-cos(u./v)./v.*exp(u+v)).^2+4.*cos(u./v).^2./v.^2.*u.^2+(v.*exp(u+v)-exp(u+v).*u).^2).^(1./2) dudv with u limits 1 to 2 and v limits 3 to 4 ans = 328.1465 % Now I will show how to control for error in your answer. % First, you can have ezint241 estimate the error in your answer by the following method. >> [a b] = ezint241(sqrt(sum(K.*K)), u, 1,2, v, 3,4) Integrate ((-exp(u+v).*cos(u./v).*u./v.^2-cos(u./v)./v.*exp(u+v)).^2+4.*cos(u./v).^2./v.^2.*u.^2+(v.*exp(u+v)-exp(u+v).*u).^2).^(1./2) dudv with u limits 1 to 2 and v limits 3 to 4 a = 328.1474 b = -8.9847e-04 % This means that the estimated answer is 328.1474 % and the estimated error in this answer is -.00089847 or about -.0009 % So probably a more accurate answer is 328.1474 - .0009 = 328.1465 % You can also specify an error tolerance in the last entry when calling ezint241. % The following asks ezint241 to have error less than .000001 times 1 plus the absolute % value of the integral. % In this case we would then be asking for an error less than about .00033 >> [a b] = ezint241(sqrt(sum(K.*K)), u, 1,2, v, 3,4,.000001) Integrate ((-exp(u+v).*cos(u./v).*u./v.^2-cos(u./v)./v.*exp(u+v)).^2+4.*cos(u./v).^2./v.^2.*u.^2+(v.*exp(u+v)-exp(u+v).*u).^2).^(1./2) dudv with u limits 1 to 2 and v limits 3 to 4 a = 328.1465 b = 1.6085e-05 % So the estimate error is .000016 and 328.1465 looks pretty good for the answer. % ezint241 uses some randomness in its method so your answer may vary if you invoke it again. >> [a b] = ezint241(sqrt(sum(K.*K)), u, 1,2, v, 3,4,.000001) Integrate ((-exp(u+v).*cos(u./v).*u./v.^2-cos(u./v)./v.*exp(u+v)).^2+4.*cos(u./v).^2./v.^2.*u.^2+(v.*exp(u+v)-exp(u+v).*u).^2).^(1./2) dudv with u limits 1 to 2 and v limits 3 to 4 a = 328.1465 b = 2.6165e-05 % Same answer above, but slightly different error.