This preview shows page 1. Sign up to view the full content.
Unformatted text preview: Homework 6 Lectures 21‐24 Math 344 Ohio University Spring Quarter 09’‐10’ Bob Benjamin Astrom Elbert Note: In the problems represented in this homework the percent error has values that are negative. The percent error is always positive however for some of the problems it is required to examine which direction (positive or negative) the error is. For this reason the absolute value has been removed from the percent error. 21.2) function [L R T] = myints(f,a,b,n) % Math 344 % Homework Problem 21.2 % Written by Astrom and Elbert % % Calculates the area under a function curve (integration) % using the Left and Right endpoint rules and the Trapezoid Rule % % inputs f, a, b and n % f is the function % a is the start point of integration % b is the end point of integration % n is number of intervals % The intervals have a start and end point which will % require n+1 points to make n intervals L = 0; % Set L to begin integration at 0 R = 0; % Set R to begin integration at 0 T = 0; % Set T to begin integration at 0 delta = (ba)/n; % Width of interval x = linspace(a,b,n+1); % Vector of n+1 points y = f(x); % Set function to variable y for i = 1:n % Loop to calculate area L=L+y(i)*delta; % Calculate using Left point R=R+y(i+1)*delta; % Calculate using Right point T=T+(y(i)+y(i+1))/2*delta; % Calculate using Trapezoid end % End calculation loop >> f = inline('sqrt(x)','x'); >> a = 1; >> b = 2; >> n = 4; >> [L R T] = myints(f,a,b,n) L = 1.166413628918445 R = 1.269967019511719 T = 1.218190324215082 >> n = 100; >> [L R T] = myints(f,a,b,n) L = 1.216879128301471 R = 1.221021263925201 T = 1.218950196113336 The percent error was calculated by adding a few extra lines of code on to the program. The extra code uses symbolic variables to calculate the definite integral with bounds “a” and “b” and then the percent error was calculated using the definite integral as the accepted value. function [L R T] = myintslong(f,a,b,n) % Math 344 % Homework Problem 21.2 % Writen by Astrom and Elbert % % Calculates the area under a function curve (integration) % using the Left and Right endpoint rules and the Trapezoid Rule % as well as the percent error associated with each method % inputs f, a, b and n % f is the function % a is the start point of integration % b is the end point of integration % n is number of intervals % The intervals have a start and end point which will % require n+1 points to make n intervals format long L = 0; % Set L to begin integration at 0 R = 0; % Set R to begin integration at 0 T = 0; % Set T to begin integration at 0 delta = (ba)/n; % Width of iterval x = linspace(a,b,n+1); % Vertor of n+1 points y = f(x); % Set function to variable y for i = 1:n % Loop to calculate area L=L+y(i)*delta; % Calculate using Left point R=R+y(i+1)*delta; % Calculate using Right point T=T+(y(i)+y(i+1))/2*delta; % Calculate using Trapezoid end % End calculation loop syms x y % Allow x and y to be symbolic true = int(f(x),x,a,b); % Definite integral true answer accepted = double(true) % Convert answer to decimal L_percent_error = double(((Ltrue)/true)*100) % Error L R_percent_error = double(((Rtrue)/true)*100) % Error R T_Percent_error = double(((Ttrue)/true)*100) % Error T >> f = inline('sqrt(x)','x'); >> a=1; >> b=2; >> n=100; >> [L R T] = myintslong(f,a,b,n) accepted = 1.218951416497460 L_percent_error = ‐ 0.170005807281802 R_percent_error = 0.169805572209664 T_Percent_error = ‐1.001175360687180e‐004 L = 1.216879128301471 R = 1.221021263925201 T = 1.218950196113336 22.2) function [M] = mymidpoint(f,a,b,n) % Math 344 % Homework Problem 22.2 % Writen by Astrom and Elbert % % Calculates the area under a function curve (integration) % using the Midpoint Rule % as well as the percent error associated with the method % % inputs f, a, b and n % f is the function % a is the start point of integration % b is the end point of integration % n is number of intervals % format long M = 0; % Set M to begin intergation at 0 delta = (ba)/n; % Width of interval halfdelta = delta/2; % Calculate half width x = linspace(a+halfdelta,bhalfdelta,n); % Vector of n points y = f(x); % Set function to variable y for i = 1:n % Loop to calculate area M = M + (y(i)*delta); % Calculate using Midpoint Rule end % End calculation loop syms x y % Allow x and y to be symbolic true = int(f(x),x,a,b); % Definite integral true answer accepted = double(true) % Convert answer to decimal M_percent_error = double(((Mtrue)/true)*100) % Error M >> f = inline('sqrt(x)','x'); >> a=1; >> b=2; >> n=4; >> [M] = mymidpoint(f,a,b,n) accepted = 1.218951416497460 M_percent_error = 0.031168549590687 M = 1.219331345974198 >> n = 100; >> [M] = mymidpoint(f,a,b,n) accepted = 1.218951416497460 M_percent_error = 5.005863613881017e‐005 M = 1.218952026687914 If we compare the percent errors calculated in 21.2 associated with the Left, Right endpoint rules and Trapezoid rules to the percent error associated with the Midpoint percent error we will see that the Midpoint rule is far more accurate than the Left and Right endpoint methods and about twice as accurate as the Trapezoid rule or ½ErrorT = ErrorM however the signs of the errors are opposite. L_percent_error = R_percent_error = T_percent_error = M_percent_error = 22.3) function [S] = mysimpson(f,a,b,n) % Math 344 % Homework Problem 22.3 % Writen by Astrom and Elbert % % Calculates the area under a function curve (integration) % using the Simpson Rule % as well as the percent error associated with the method % % inputs f, a, b and n % f is the function % a is the start point of integration % b is the end point of integration % n is number of intervals % if rem(n,2) ~= 0 % Check to see that n is even error('n must be even') % Error message end S = 0; % Set S to begin integration at 0 x = linspace(a,b,n+1); % Vector of n+1 points y = f(x); % Set function to variable y w = ones(n+1,1); % Fill w with ones for i = 2:n % Loop of inner (nonendpoints) n's if rem(i,2)==0 % Check to see if i is even w(i)=4; % If n is even weight equals 4 else % If i is odd w(i)=2; % weight equals 2 end % End of loop S = S + w(i)*y(i); % Calculates sum of weighted points end % End loop S = (S + y(1) + y(n+1))*((ba)/(3*n)); % Add enpoints and multiply by % interval width divided by 3n. syms x y % Allow x and y to be symbolic true = int(f(x),x,a,b); % Definite integral true answer accepted = double(true) % Convert answer to decimal S_percent_error = double(((Strue)/true)*100) % Error S ‐ 0.170005807281802 0.169805572209664 ‐1.001175360687180e‐004 5.005863613881017e‐005 >> f = inline('sqrt(x)','x'); >> a=1; >> b=2; >> n=4; >> [S] = mysimpson(f,a,b,n) accepted = 1.218951416497460 S_percent_error = ‐5.135266499708094e‐004 S = 1.218945156857086 >> n=100; >> [S] = mysimpson(f,a,b,n) accepted = 1.218951416497460 S_percent_error = ‐1.406795518215536e‐009 S = 1.218951416480312 When the percent error of the Simpson method is compared to the previous Endpoint, Trapezoid and Midpoint methods the Simpson method is found to be far more accurate. 23.1) Since 3 – (‐3) = 6, selecting m = 60 would let h = .1 and because 2 – (‐2) = 4, letting n = 40 would also let k = .1; thereby producing the following the following: >> [X Y] = meshgrid(‐3:.1:3, ‐2:.1:2) >> f = inline('x.*exp(‐x.^2 ‐ y.^2)', 'x', 'y') >> Z = f(X, Y) >>mesh(X, Y, Z) 23.2) % mytriangles % % Math 344 % Homework Problem 23.2 % Writen by Astrom and Elbert % % Program to produce a triangulation. % V contains vertices, which are (x,y) pairs V = [ 2 2 ; 2 2 ; 4 0 ; 4 0 ; 2 2 ; 2 2 ; 0 1 ; 0 1 ; 1 0 ; 1 0 ] % x, y are row vectors containing coordinates of vertices x = V(:,1) y = V(:,2) T = delaunay(x,y) % Assign the triangles trimesh(T,x,y) % Plot the output >> mytriangles V = ‐2 2 2 2 ‐4 0 4 0 2 ‐2 ‐2 ‐2 0 1 0 ‐1 ‐1 0 1 0 x = ‐2 2 ‐4 4 2 ‐2 0 0 ‐1 1 y = 2 2 0 0 ‐2 ‐2 1 ‐1 0 0 T = 9 1 3 5 8 6 6 9 3 8 9 6 10 5 4 4 2 10 10 8 5 7 2 1 1 9 7 7 10 2 8 10 7 7 9 8 The following is the plot created with these coordinates: Figure 23.2 A) 2 dimensional plan view of convex irregular region The following is the program that creates the triangles, but also creates a z‐dimension so that a 3‐ dimensional plot is created (this could possibly represent a container or a scoop of some sorts): % mytrianglesplot % % Math 344 % Homework Problem 23.2 % Writen by Astrom and Elbert % % Program to produce a triangulation. % V contains vertices, which are (x,y) pairs V = [ 2 2 ; 2 2 ; 4 0 ; 4 0 ; 2 2 ; 2 2 ; 0 1 ; 0 1 ; 1 0 ; 1 0 ] % x, y are row vectors containing coordinates of vertices x = V(:,1) y = V(:,2) f = inline('abs(sin(x.*y)).^(3/2)','x','y'); % Function to calculate z z = f(x,y) % Calculation of z T = delaunay(x,y) % Assign the triangles trimesh(T,x,y,z) % Plot the output >> mytrianglesplot V = ‐2 2 2 2 ‐4 0 4 0 2 ‐2 ‐2 ‐2 0 1 0 ‐1 ‐1 0 1 0 x = ‐2 2 ‐4 The following is the plot created by mytrianglesplot: 4 2 ‐2 0 0 ‐1 1 y = 2 2 0 0 ‐2 ‐2 1 ‐1 0 0 z = 0.65837576047485 0 0.65837576047485 0 0 0 0.65837576047485 0 0.65837576047485 0 0 0 0 0 T = 9 1 3 5 8 6 6 9 3 8 9 6 10 5 4 4 2 10 10 8 5 7 2 1 1 9 7 7 10 2 8 10 7 7 9 8 Figure 23.2 B) Isometric view of convex irregular region with z values 24.1) The following is mylowerleft with a few extra lines of code to calculate the percent error: function I = mylowerleft(f,a,b,c,d,m,n) % % Math 344 % Homework Problem 23.2 % Writen by Astrom and Elbert % % Computes an approximate integral of a function of two variables f(x,y) % Uses the lowerleft corner to evaluate. % Inputs: % a,b  define the interval in x, namely a<x<b % c,d  define the interval in y, namely c<x<d % m  the number of (evenly spaced) intervals to use in x % n  the number of (evenly spaced) intervals to use in y % Output: the approximate value of the integral % h = (ba)/m; % step size in x direction k = (dc)/n; % step size in y direction [X,Y] = meshgrid(a:h:b,c:k:d); % sets up a grid Z=f(X,Y); mesh(X,Y,Z) % plot, just for fun I = 0; % start adding at zero for i = 1:m % x loop for j = 1:n % y loop I = I + Z(j,i); % add up values of f(x,y) % note that x corresponds to the column index, so i is second end end I = I*h*k; answer_true_integral = dblquad(f,a,b,c,d) percent_error=(I  answer_true_integral)/answer_true_integral*100 >> a = 1; >> b = 5; >> c = 0; >> d = 2; >> m = 10; >> n = 18; >> I = mylowerleft(f,a,b,c,d,m,n) answer_true_integral = 12.797482353249382 percent_error = ‐8.083879029107308 I = 11.762949361041347 Figure 24.1 A) Isometric view of Lower Left Figure 24.1 B) Plan view of Lower Left The following is the program modified to calculate integral using the center point method: function I = mycenterpoint(f,a,b,c,d,m,n) % % Math 344 % Homework Problem 24.1 % Writen by Astrom and Elbert % % Computes an approximate integral of a function of two variables f(x,y) % Uses the center point to evaluate. % Inputs: % a,b – define the interval in x, namely a<x<b % c,d – define the interval in y, namely c<x<d % m  the number of (evenly spaced) intervals to use in x % n  the number of (evenly spaced) intervals to use in y % Output: the approximate value of the integral h = (ba)/m; % step size in x direction k = (dc)/n; % step size in y direction [X,Y] = meshgrid((a + h/2):h b – h/2),(c + k/2):k d – k/2)); % sets up a grid Z=f(X,Y); mesh(X,Y,Z) % plot, just for fun I = 0; % start adding at zero for I = 1:m % x loop for j = 1:n % y loop I = I + Z(j,i); % add up values of f(x,y) % note that x corresponds to the column index, so I is second end end I = I*h*k; answer_true_integral = dblquad(f,a,b,c,d) % Double integration percent_error=(I – answer_true_integral)/answer_true_integral*100 % percent % error >> f = inline('sqrt(x.*y)','x','y'); >> a = 1; >> b = 5; >> c = 0; >> d = 2; >> m = 10; >> n = 18; >> I = mycenterpoint(f,a,b,c,d,m,n) answer_true_integral = 12.797482353249382 percent_error = 0.137017495840315 I = 12.815017143100411 Figure 24.1 C) Isometric View of the graph of mycenterpoint Figure 24.1 D) Plan view of the graph of mycenterpoint Based on the percent error of the two methods it can be concluded that the Center Point Method is far more accurate than the Lower Left Point Method. 24.2) function I = mydblsimp(f,a,b,c,d,m,n) % % Math 344 % Homework Problem 24.2 % Writen by Astrom and Elbert % % Produces the m by n matrix of weights for Simpson's rule % for double integrals % Inputs: f  function to be approximated % a  lower bound of x domain % b  upper bound for x domain % c  lower bound for y domain % d  upper bound for y domain % m  number of intervals in the row direction must be even % n  number of intervals in the column direction must be even % Output: W  a (m+1)x(n+1) matrix of the weights if rem(m,2)~=0  rem(n,2)~=0 % Check that matrix is even error('m and n must be even') % Error message end u = mysimpweights(m); % X weight set up v = mysimpweights(n); % Y weight set up W = u*v'; % Creates weighted matrix h = (ba)/m; % step size in x direction k = (dc)/n; % step size in y direction [X,Y] = meshgrid(a:h:b,c:k:d); % sets up a grid Z=f(X,Y); mesh(X,Y,Z) % plot, just for fun I = 0; % start adding at zero for i = 1:m+1 % x loop for j = 1:n+1 % y loop I = I + W(i,j).*Z(j,i); % add up values of f(x,y) % multiplied by the corresponding % Simpson coefficient note that x % corresponds to the % column index, so i is second end end I = (1/9)*I*h*k; answer_true_integral = dblquad(f,a,b,c,d) % Double integration percent_error=(I  answer_true_integral)/answer_true_integral*100 % percent % error >> f = inline('sqrt(x.*y)','x','y'); >> a = 1; >> b = 5; >> c = 0; >> d = 2; >> m = 10; >> n = 18; >> I = mydblsimp(f,a,b,c,d,m,n) answer_true_integral = 12.797482353249382 Percent_error = ‐0.160082562014088 I = 12.776995815625000 Figure 24.2 A) Isometric view of the graph for mydblsimp Figure 24.2 B) Plan view of the graph for mydblsimp The percent error for the Simpson method is slightly higher than the Center Point method. The Simpson method returns a percent error of 0.16% where as the center point method has an error of 0.13%. In this case the Center Point method is more accurate than the Simpson method 24.3) The following program, mytrapmethod (code listed second), uses a function program called mytrapweights, whose code is: function w = mytrapweights(n) % % Math 344 % Homework Problem 24.3 part 1 % Writen by Astrom and Elbert % % computes the weights for the trapezoid rule % Input: n  the number of intervals, must be even % Output: a vector with the weights, length n+1 if rem(n,2) ~= 0 % Check for odd numbers error('n must be even') % Error message end w = ones(n+1,1); % Enter ones for i = 2:n % Loop for nonend points w(i)=2; % Weights for inner (nonend point) points end >> f = inline('sqrt(x.*y)','x','y'); >> a = 1; >> b = 5; >> c = 0; >> d = 2; >> m = 10; >> n = 18; >> w = mytrapweights(n) w = 1 2 2 2 function I = mytrapmethod(f,a,b,c,d,m,n) % % Math 344 % Homework Problem 24.3 part 2 % Writen by Astrom and Elbert % % Produces the m by n matrix of weights for Trapezoid rule 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 % for double integrals % Inputs: f  function to be approximated % a  lower bound of x domain % b  upper bound for x domain % c  lower bound for y domain % d  upper bound for y domain % m  number of intervals in the row direction. % must be even. % n  number of intervals in the column direction. % must be even. % Output: I  finite approximation of the trapezoid method for integrals if rem(m,2)~=0  rem(n,2)~=0 error('m and n must be even') end u = mytrapweights(m); v = mytrapweights(n); W = u*v' h = (ba)/m; k = (dc)/n; [X,Y] = meshgrid(a:h:b,c:k:d); Z=f(X,Y); I = 0; for i = 1:m+1 for j = 1:n+1 I = I + W(i,j).*Z(j,i); % step size in x direction % step size in y direction % sets up a grid % start adding at zero % x loop % y loop % add up values of f(x,y) % multiplied by the corresponding % Trapezoid coefficient % note that x corresponds to the % column index, so i is second end end I = (1/4)*I*h*k; >> f = inline('sqrt(x.*y)','x','y'); >> a = 1; >> b = 5; >> c = 0; >> d = 2; >> m = 10; >> n = 18; >> I = mytrapmethod(f,a,b,c,d,m,n) W = Columns 1 through 12 1 2 2 2 2 2 2 2 2 2 2 2 2 4 4 4 4 4 4 4 4 4 4 4 2 4 4 4 4 4 4 4 4 4 4 4 2 4 4 4 4 4 4 4 4 4 4 4 2 4 4 4 4 4 4 4 4 4 4 4 2 4 4 4 4 4 4 4 4 4 4 4 2 4 4 4 4 4 4 4 4 4 4 4 2 4 4 4 4 4 4 4 4 4 4 4 2 4 4 4 4 4 4 4 4 4 4 4 2 4 4 4 4 4 4 4 4 4 4 4 1 2 2 2 2 2 2 2 2 2 2 2 Columns 13 through 19 2 2 2 2 2 2 1 4 4 4 4 4 4 2 4 4 4 4 4 4 2 4 4 4 4 4 4 2 4 4 4 4 4 4 2 4 4 4 4 4 4 2 4 4 4 4 4 4 2 4 4 4 4 4 4 2 4 4 4 4 4 4 2 4 4 4 4 4 4 2 2 2 2 2 2 2 1 I = 12.740803502534769 ...
View
Full
Document
This note was uploaded on 01/15/2011 for the course MATH 345 taught by Professor Staff during the Spring '08 term at Ohio State.
 Spring '08
 Staff

Click to edit the document details