HW 8 Solution - Z(:,3) = t(n-3:n).^2; Z(:,4) = t(n-3:n).^3;...

Info iconThis preview shows pages 1–5. Sign up to view the full content.

View Full Document Right Arrow Icon
HW 8 Solution Part 1. A-C. load cvt.dat t = cvt(:,1); c = cvt(:,2); % Compute mol/min, c is in mol/L F = 1000; % L/min dMdt = F*c; % mol/min %plot(t,dMdt) % I.A. Trapezoid rule is same as NC p = 1 n = length(t); dM1 = zeros(n,1); M1 = zeros(n,1); for i = 1:n-1 dM1(i) = (t(i+1)-t(i))*(dMdt(i)+dMdt(i+1))/2; M1(i+1) = M1(i) + dM1(i); end % I.B. Newton-Cotes p = 2 dM2 = zeros(n,1); M2 = zeros(n,1); Z = zeros(3); for i = 1:2:n-4 y = dMdt(i:i+2); Z(:,1) = ones(3,1); Z(:,2) = t(i:i+2); Z(:,3) = t(i:i+2).^2; a = inv(Z'*Z)*Z'*y; dM2(i) = a(1)*(t(i+2)-t(i)) + a(2)/2*(t(i+2)^2-t(i)^2) + . .. a(3)/3*(t(i+2)^3 - t(i)^3); M2(i+2) = M2(i) + dM2(i); end % Finish off final piece it is exists with p = 3 if (mod(n,2) == 0) y = dMdt(n-3:n); Z = zeros(4); Z(:,1) = ones(4,1); Z(:,2) = t(n-3:n);
Background image of page 1

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon
Background image of page 2
Background image of page 3

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon
Background image of page 4
Background image of page 5
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: Z(:,3) = t(n-3:n).^2; Z(:,4) = t(n-3:n).^3; a = inv(Z'*Z)*Z'*y; dM2(n-3) = a(1)*(t(n)-t(n-3)) + a(2)/2*(t(n)^2-t(n-3)^2) + . .. a(3)/3*(t(n)^3 - t(n-3)^3) + a(4)/4*(t(n)^4 - t(n-3)^4); M2(n) = M2(n-3) + dM2(n-3); end % I.C. Plot figure(2) plot(t,M1,t,M2,'x') xlabel('time (min)') ylabel('total product (mol)') Pat 2. A-C. % Part II: Differentiate % II.A. D1 = zeros(n,1); for i = 1:n-1 D1(i) = (M1(i+1)-M1(i))/(t(i+1)-t(i)); end D1(n) = (M1(n)-M1(n-1))/(t(n)-t(n-1)); % II.B. D2 = zeros(n,1); for i = 2:n-1 A = [t(i+1)-t(i) 0.5*(t(i+1)-t(i))^2; t(i-1)-t(i) 0.5*(t(i-1)-t(i))^2]; b = [M1(i+1)-M1(i); M1(i-1)-M1(i)]; x = inv(A)*b; D2(i) = x(1); end figure(1) plot(t,c,t,D1/F,'--',t,D2/F,'o') xlabel('time (min)') ylabel('concentration (mol/L)')...
View Full Document

This note was uploaded on 05/10/2008 for the course CHBE 2120 taught by Professor Gallivan during the Summer '07 term at Georgia Institute of Technology.

Page1 / 5

HW 8 Solution - Z(:,3) = t(n-3:n).^2; Z(:,4) = t(n-3:n).^3;...

This preview shows document pages 1 - 5. Sign up to view the full document.

View Full Document Right Arrow Icon
Ask a homework question - tutors are online