MATLAB

SA 11
num_int_ex_moodle.m

% numerical integration example using trap rule. function [Da,Db] = num_int_ex_moodle() tic % givens t1 = 0; tf = 50; %% calculate distance using num int for different values of N N = [2^4; 2^8; 2^12; 2^16; 2^20]; % differet values of n for k=1:length(N) % loop over each value of n, will create a vector for distance for each value n = N(k); % the particular value of n for this iteration t = linspace(t1,tf,n+1)'; % equally spaced vector of time between t1 and tf h = (tf-t1)/n; % segment size % solve using a for loop, Da is eatimate of distance for this approach tic % keep track of time Da(k,1) = vel(t1); % v(t) at first time step, note Da(k,1) is the estimate of D for the kth value of n for j=2:n % loop over points 2-n and add to the integral Da(k,1) = Da(k,1) + 2*vel(t(j)); end Da(k,1) = Da(k,1) + vel(tf); % add final value of v(t) at last time step Da(k,1) = Da(k,1)*h/2; % now scale sum by h/2 for trapezoid method tma(k)=toc; % finish of timing % solve with a vectorized approach not a for loop, Db is eatimate of distance for this approach tic Db(k,1) = vel(t1); % v(t) at first time step Db(k,1) = Db(k,1) + 2*sum(vel(t(2:n,1))); % sum all values of v(t) for 2-n without loop Db(k,1) = Db(k,1) + vel(tf); % add final value of v(t) at last time step Db(k,1) = Db(k,1)*h/2; % now scale sum by h/2 for trapezoid method tmb(k)=toc; % finish of timing if k>1 % calculate approximate error for increasing values of n e(k-1) = 100*abs((Db(k,1)-Db(k-1,1))/Db(k,1)); end end % plot the estimate of the integral (note Da=Db so only need one of them) figure(2) plot(N,Db) axis([0 2^21 0 ceil(Db(k)/10^3)*10^3]) ylabel('Distance') xlabel('N') % plot approx error on a loglog plot figure(3) loglog(N(2:length(N)),e) ylabel('Approximate Error') xlabel('N') % plot timimg on a loglog plot figure(4) loglog(N,tma) hold on loglog(N,tmb,'r') legend('trap, forloop','trap, vec') ylabel('Comp time') xlabel('N') toc end % function that calculate V as a function of t function [V] = vel(t) g=9.81; m=70; c = 12.5; V = g*m/c*(1-exp(-c/m.*t)); end