Numerical Methods with Chemical Engineering Applications 1st Edition

Numerical methods with chemical engineering applications Solutions Manual Kevin D. Dorfman and Prodromos Daoutidis Contents List of illustrations page iv List of tables v 1 Introduction 1 2 Linear equations 43 3 Nonlinear equations 115 4 Initial value problems 208 5 Dynamical systems 276 6 Boundary value problems 346 7 Partial differential equations 402 8 Interpolation and integration 468 Illustrations 2.1 Process schematic for Problem 2.80 61 5.1

Introduction 25 15 for l = 1:n 16 A(k,l) = (k+l)/(k*l); 17 end 18 end 19 At = transp(A,m,n);20 C = mult(A,At,m,n);21 tr = tra(C,m);22 row entry = (i-1)*m max + j;23 output(row entry,1) = m;24 output(row entry,2) = n;25 output(row entry,3) = tr;26 end27 end28 dlmwrite('s11c1p1 output.txt',output); 19 At = transp(A,m,n);20 C = mult(A,At,m,n);21 tr

Linear equations 65 (2.8) Cofactor expansion for i=3 has the most zero entries, giving det A=(1)(1)(3+1) det 2 1 0 1 3214 7234 3 4 1 2 3 4 1 2 3 4 1 2 +(0)(1)(3+2) det 8 2 3 44412+(1)(1)(3+3) det 3 2 0 153148734 +(0)(1)(3+2) det 8 2 3 44412+(1)(1)(3+3) det 3 2 0 153148734 +(0)(1)(3+2) det 8

Linear equations 85 We then solve241055x1x2=00 We then solve241055x1x2=00 We then solve 241 055 x1 x2 = 0 0 004 x3 1 which gives x1x2=3/81/4 x1x2=3/81/4 x1 x2 = 3/8 1/4 x3 1/4 So the inverse is A1=9/20 1/20 1/4 A1=9/20 1/20 1/4 A1= 9/20 1/20 1/4 11/40 1/40 3/8 11/40 1/40 3/8 11/40 1/40 3/8 1/4 1/4 1/4 or, if you prefer nicer numbers A1=140 18 21010 10 10 A1=140 18 21010 10 10 A1=1 40 18 210 10 10 10 11 1 15 11 1 15 11 1 15 (2.43)

Linear equations 105 5A = [6, 1; 4, -8]; 6b = [13;0]; 7x solve = A\b; 8 9x = [5,0];10 [xplot,yplot,err plot,k] = jacobi(A,x,b);11 xplot = xplot;12 yplot = yplot;1314 %make the original figure15 h = figure;16 plot(xplot,yplot,'-ob')17 xlabel('$x$','FontSize',14), ylabel('$y$','FontSize',14)18 title('Solution to s12c4p2','FontSize',14) 9x = [5,0];10 [xplot,yplot,err plot,k] = jacobi(A,x,b);11 xplot = xplot;12 yplot = yplot;1314

Nonlinear equations 139 The initial guess x(0) = 0.2 0.2 0.2 0.2 0.2 converges to x()= 1 0 1 in 6 iterations.The Matlab program is:1function s10c4p32clc3x zero = [0.2;0.2;0.2];4%check to see if input is in column form5n rows = size(x zero,1);6if n rows == 17x zero = x zero.';8end in 6 iterations.The Matlab program is:1function s10c4p32clc3x zero = [0.2;0.2;0.2];4%check to see if input is

Nonlinear equations 159 65 fprintf('Temperature difference at inlet is ... \t%8.6f\n',T1-T1p) 66 fprintf('Temperature difference at outlet is ... \t%8.6f\n\n',T2-T2p) 67 Qinner = w*Cp*(T1-T2);68 Qouter = wp*Cpp*(T1p-T2p);69 fprintf('The Q value of the inner fluid is ...\t%8.6f\n',Qinner)70 fprintf('The Q value of the outer fluild is ...\t%8.6f\n',Qouter)71 Tlm = ((T1-T1p)-(T2-T2p))/log((T1-T1p)/(T2-T2p));72 Qtransfer = U*A*Tlm;73 fprintf('The Q value from heat

Nonlinear equations 179 53 xlabel('$x$','FontSize',14) 54 ylabel('$y$','FontSize',14) 55 title('Solution to s13c6p3','FontSize',14) 56 saveas(g,'s13c6p123 solution figure3.eps') 5758 function out = getR(P1sat,P2sat,A,x,z)59 y = z(1); P = z(2); %unpack the vector60 R(1,1) = x*P1sat*exp(A*(1-x)2) - y*P;61 R(2,1) = (1-x)*P2sat*exp(A*x2) - (1-y)*P;62 out = R;6364 function out = getJ(z)65 y = z(1); P = z(2); %unpack the vector66

Nonlinear equations 199 325 %this is a vapor 326 P = vdw(V,T,a,b); 327 else 328 %this is on the horizontal part 329 P = Psat;330 end331 P low(i) = P;332 end333334 h=figure;335 plot(Vplot,Pvdw,'-or',Vplot,P low,'-ob')336 legend('EOS values','Isotherm')337 axis([Vmin,Vmax,Pmin,Pmax])338 xlabel('$\underline{V}$ (L/mol)','FontSize',14) 329 P = Psat;330 end331 P low(i) = P;332 end333334 h=figure;335 plot(Vplot,Pvdw,'-or',Vplot,P low,'-ob')336 legend('EOS values','Isotherm')337 axis([Vmin,Vmax,Pmin,Pmax])338 xlabel('$\underline{V}$

Initial value problems 231 For the third equation, let y3=xso we get d dxy2sin y3y24 dxy2sin y3y24 dx y2 sin y3 y2 4 =,0= =,0= = , 0 = y0 y1 y31y30 y31y30 y3 1 y3 0 y1 y2 y0 y1 1 0 For the last equation, we use the product rule to gety =1(y)2yConverting to a system of equations yieldsddx"y0y1#="y11y2y1#0="21# For the last equation, we use the product rule to gety =1(y)2yConverting to a system of equations yieldsddx"y0y1#="y11y2y1#0="21# For the last

Initial value problems 251 51 function out = feval(x,y) 52 out = x/sqrt(y)+sin(pi*x*y);The output file is:1 1.5 2 2.5 3 3.5 4 4.5 512345678xySolut ion for s 12c 8p2h= 0. 01h=hh=h+ 0. 1(4.36) The files for this problem are contained in the folder s11c5p1 matlab.The Matlab script is:1function s11c5p12close all3clc4set(0,'defaulttextinterpreter','latex')56ncalcs = 100;78n

Initial value problems 271 of equations is dAdt =k1AB +k2CdBdt =k1AB +k2Ck3BC +k4DEdCdt =k1AB k2Ck3BC +k4DEdDdt = +k3BC k4DEdEdt = +k3BC k4DE k5EdF dAdt =k1AB +k2CdBdt =k1AB +k2Ck3BC +k4DEdCdt =k1AB k2Ck3BC +k4DEdDdt = +k3BC k4DEdEdt = +k3BC k4DE k5EdF dA dt =k1AB +k2C dB dt =k1AB +k2Ck3BC +k4DE dC dt =k1AB k2Ck3BC +k4DE dD dt = +k3BC k4DE dE dt = +k3BC k4DE

Dynamical systems 305 2clc 3close all 4set(0,'defaulttextinterpreter','latex') 5 6h = figure;7hold on89x = [0.5,0];10 makeplot(x,'-b')1112 x = [-2,3];13 makeplot(x,'-r')1415 hold off 6h = figure;7hold on89x = [0.5,0];10 makeplot(x,'-b')1112 x = [-2,3];13 makeplot(x,'-r')1415 hold off 6h = figure; 7hold on 8 9x = [0.5,0]; 10 makeplot(x,'-b') 11 12 x = [-2,3]; 13 makeplot(x,'-r') 14 15 hold off 16 17 xlabel('$x 1$','FontSize',14) 18 ylabel('$x 2$','FontSize',14) 19 20 saveas(h,'s12c9p3 solution figure.eps','psc2') 212223 function makeplot(x,plotformat)24 h

Dynamical systems 325 72 function y = RK4(y0, h, Nt) 737475 y = zeros(Nt+1, 2);76 y(1, :) = y0;7778 for i = 1:Nt7980 yi = y(i, :);8182 dy1 = h*f( yi );83 dy2 = h*f( yi + 0.5*dy1 );84 dy3 = h*f( yi + 0.5*dy2 );85 dy4 = h*f( yi + dy3

Dynamical systems 3450.1 0.08 0.06 0.04 0.02 0 0.02 0.04 0.06 0.08 0.10.10.080.060.040.0200.020.040.060.080.1xvS ol ut ion t o s11c 7p 1 Dynamical systems 3450.1 0.08 0.06 0.04 0.02 0 0.02 0.04 0.06 0.08 0.10.10.080.060.040.0200.020.040.060.080.1xvS ol ut ion t o s11c 7p 1 Dynamical systems 345 0.1 0.08 0.06 0.04 0.02 0 0.02 0.04

372 Boundary value problems 57 A(1,1) = 1; b(1) = 1; 58 59 %right BC 60 A(n,n) = 1; b(n) = 0; 6162 %interior nodes63 for i = 2:n-164 A(i,i-1) = 1;65 A(i,i) = -2 - dx2;66 A(i,i+1) = 1;67 end6869 %solve equation70 A = sparse(A);71 c=A\b;Note that we used sparse matrices to force Matlab

392 Boundary value problems 37 k2 = k2 + dk; 38 else 39 fprintf('Amazingly got r = 1 to numerical ... precision.\n') 40 end4142 err = abs(r-1);43 end44 fprintf('\nConverged to within tolerance of %2.1e.\n',tol)45 fprintf('\n****************************************\n')464748 function out = compute c(Da,Db,k1,k2)49 %setup of the grid 40 end4142 err = abs(r-1);43 end44 fprintf('\nConverged to within tolerance of %2.1e.\n',tol)45

Partial differential equations 423 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 P os i ti on , x C on c e nt r at i on , c(x , t ) S ol u ti o n t o s 13 c 11p 3 0. 0005 0. 005 0. 02 0. 05 0. 1 0.

Partial differential equations 443 ck+n. The new equations for y=0 are 2ck+n+ck+14ck+ck1=0 2ck+n+ck+14ck+ck1=0 2ck+n+ck+14ck+ck1=0 We also need to change the lower corners to no-ux boundary conditions. This leads to cn+1+c22c1=0andc2n+cn12cn=0 cn+1+c22c1=0andc2n+cn12cn=0 cn+1+c22c1=0 and c2n+cn12cn=0 The no-ux corner boundary condition is discussed in the last recitation. The exact solution is c=1. The Matlab file for this problem is 1function s15h12p2 2clc 3close all 4set(0,'defaulttextinterpreter','latex') 5 6%setup grid 7n = 101; 8dx = 1/(n-1); 9A

Partial differential equations 463 and the ratio is Q2=2Ldx Q2=2Ldx Q2 =2L dx Qr 2 TrZ1 1/2 d2 y!y=L The Matlab script for this problem is: 1function s13c12p3 2clc 3close all 4set(0,'defaulttextinterpreter','latex')56%set the problem parameters7k = 10; %ratio of conductivities k1/k28L = 0.5; %aspect ratio of the channel 4set(0,'defaulttextinterpreter','latex')56%set the problem parameters7k = 10; %ratio of conductivities k1/k28L = 0.5; %aspect ratio of the channel 4set(0,'defaulttextinterpreter','latex') 5 6%set the problem

488 Interpolation and integration 488 Interpolation and integration 488 Interpolation and integration (8.11) The files for this problem are contained in the folder s12c7p2 matlab. The Matlab script is: 1function s12c7p2 2clc 3close all 4set(0,'defaulttextinterpreter','latex')56%make a plot of the function for my own check7fplot(@feval,[-1,1]);8xlabel('$x$','FontSize',16), ylabel('$f(x)$','FontSize',16)910 %compute each of the numerical integrals1112 %trapezoidal rule13 f = get points(2); 4set(0,'defaulttextinterpreter','latex')56%make a