Nonlinear equations 139
The initial guess
x(0) =
0.2
0.2
converges to
x(∗)=
1
0
1
9
10 x = x zero;
11 err = 1000;
12 i = 0;
13 while err >0.001
24
25 function out = res(w)
26 x = w(1); y = w(2); z = w(3);
27 R(1) = (xˆ3-xˆ2)/z + sin(pi*y);
28 R(2) = x*yˆ3 – cos(pi*z)-z;
39 J(2,3) = pi*sin(pi*z)-1;
40 J(3,1) = 1/z – pi*sin(pi*(x-z));
140 Nonlinear equations
41 J(3,2) = -2*exp(y);
42 J(3,3) = -x/zˆ2+pi*sin(pi*(x-z));
(3.23) The files for this problem are contained in the folder s14c6p2 matlab.
From the first equation, we either have x=0 or x=y. For x=0, the second
equation gives us y=1, which is a real root
x∗=“0
1#
1.1#
and got the following output:
0 1.4142e-01
2 3.0411e-06
4 3.6351e-23
which is quadratic convergence .
The Matlab program is:
Nonlinear equations 141
1function s14c6p2
12 fprintf(‘%2d \t %6.4e \n’,k,ek)
13
14 while ek >1e-15
15 %do Newton-Raphson
16 R = getR(x); %get residual
27 if k>99
28 fprintf(‘Did not converge\n’)
29 break
30 end
31 end
(3.24) The files for this problem are contained in the folder s14c5p12 matlab
(a) You need to first rewrite the problem the form f(x)=0,
c=7
4
Dp
and, in the Matlab program, I defined x=ǫ. So the problem I was trying to solve
is the root of
For Newton’s method, we also need
The Matlab program is
1function s14c5p1
2clc
3
4eta = 1; %cp
5vo = 0.1; %m/s
6Dp = 5; % cm
7rho = 2; %g/cmˆ3
8dP = 416; % Pa
9L = 1.5; %m
20
21 fprintf(‘Starting Newton Method\n’)
22 fprintf(‘==================================================\n’)
23 fprintf(‘k \t epsilon \t f \n’)
24 %do Newton’s method to get epsilon = x
34 printf(‘Did not converge!\n’)
35 out = -1; %dummy value
36 break
37 end
38 f = FunEval(x,a,b,c); %for while loop and next iteration
49 term2 = (-xˆ3 – 3*xˆ2*(1-x))/xˆ6;
50 out = -b*term1 – c*term2;
To see quadratic convergence, I chose an initial guess ǫ(0) =0.1 and a tolerance
of 10−12. The output to screen is
1Starting Newton Method
12 8 0.701500 -3.2950e+02
13 9 0.772616 -6.8457e+01
14 10 0.796499 -5.0679e+00
15 11 0.798564 -3.3613e-02
16 12 0.798578 -1.5037e-06
However, if you know the root, you can go back and check for quadratic conver-
gence. So a better program would be the following:
144 Nonlinear equations
1function s14c5p1 quadratic
2clc
13 Dp = Dp/100; %in m
14 eta = eta*1e-3; %Pa*s
15
16 %combine into the lumped variables
17 a = dP/L;
28 while abs(f) >1e-12
29 fprintf(‘%2d \t %8.6f \t %6.4e \n’,count,x,f)
30 df = DerEval(x,b,c);
31 x = x – f/df;
32 count = count + 1;
43
44 fprintf(‘Checking for quadratic convergence:\n’)
45 fprintf(‘==================================================\n’)
46 fprintf(‘k \t epsilon \t|x-x*| \n’)
47 x star = x; %set as the converged result for the root
57 printf(‘Did not converge!\n’)
58 out = -1; %dummy value
59 break
60 end
61 f = FunEval(x,a,b,c); %for while loop and next iteration
40 0.100000 6.9858e-01
62 0.173974 6.2460e-01
84 0.296731 5.0185e-01
95 0.382114 4.1646e-01
10 6 0.483867 3.1471e-01
16 12 0.798578 6.2108e-10
(b) I modified the program from part a so that it loops over different initial guesses
and outputs the number of iterations. The old program has become a subfunction
1function s14c5p2
2clc
3close all
4set(0,‘defaulttextinterpreter’,‘latex’)
5
15 Dp = Dp/100; %in m
16 eta = eta*1e-3; %Pa*s
17
18 %combine into the lumped variables
19 a = dP/L;
30
31 %make the plot
32 h = figure;
33 plot(x,n,‘-ob’,‘MarkerSize’,8)
34 xlabel(‘Initial Guess for $\epsilon$’,‘FontSize’,14)
45 tol = 5.844e-14;
46 fprintf(‘Starting Newton Method\n’)
47 fprintf(‘==================================================\n’)
48 fprintf(‘k \t epsilon \t f \n’)
49 %do Newton’s method to get epsilon = x
60 break
61 end
62 f = FunEval(x,a,b,c); %for while loop and next iteration
63 end
64 fprintf(‘%2d \t %8.6f \t %6.4e \n’,count,x,f)
Nonlinear equations 147
66 fprintf(‘The converge value is epsilon = %16.14f\n’,x)
67 out = count;
76 out = -b*term1 – c*term2;
You will notice that I had to set the tolerance to 5.884e-14 to get this to converge
architecture of the computer and how you write the program.
The output is
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
2
4
6
8
10
12
14
16
18
20
22
I n it ia l Gu e s s f o r ²
N umb e r o f I t e r at i on s
S ol ut io n t o s 14c 5 p2
(3.25) The files for this problem are contained in the folder s11c4p1 matlab.
The Matlab script is:
1function s11c4p1
2clc
3set(0,‘defaulttextinterpreter’,‘latex’)
13 x zero = 1;
14 tol = 1e-6;
15
16
17 %analytical derivative
28 fprintf(‘The analytical derivative converged to %8.6f after …
%2d iterations.\n’,x,iter)
29 dx output(1) = 0;
30 iter output(1) = iter;
31 x output(1) = x;
42 x = x – f/fderN(x,dx);
43 f = feval(x);
44 err = abs(f);
45 if iter >250
46 fprintf(‘Did not converge!\n’)
56 dx output(i+1) = dx;
57 iter output(i+1) = 0;
58 x output(i+1) = 0;
59 end
60 end
70 output save(:,2) = x output;
71 output save(:,3) = iter output;
72 dlmwrite(‘s11c4p1 output.txt’,output save)
73
74 function out = feval(x)
30.02,1.0809,4
40.03,1.0809,4
50.04,1.0809,4
60.05,1.0809,4
70.06,1.0809,5
90.08,1.0809,6
11 0.1,1.0809,7
13 0.12,1.0809,10
15 0.14,1.0809,15
16 0.15,1.0809,19
150 Nonlinear equations
35 0.34,1.2312,46
36 0.35,0,0
37 0.36,1.2312,223
38 0.37,1.2312,25
39 0.38,0,0
40 0.39,0,0
The file output is:
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
0
50
100
150
200
250
dx
k
S ol ut io n t o s1 1c 4p 1
Note that the case dx =0.19 converges to 1.2312 after 2939 iterations.
(3.26) The files for this problem are contained in the folder s14c6p3 matlab.
The chemical equilibrium equations we need to solve are
where EB subscripts stand for ethylbenzene, T subscripts stand for toluene, and pjis
the partial pressure of species j. If we write the unknown vector as
Nonlinear equations 151
then the Jacobian is
Psat
EB xEB
dPsat
EB
dT
The Matlab program for this problem is
1function s14c6p3
2clc
3
4x = [0.5; 75]; %initial guess for [x,T]
5R = getR(x); %get now for checking while loop
6err = norm(R);
17 delta = -J\R;
18 x = x + delta;
19 k=k+1;
20
21 %get R for next step and display error
152 Nonlinear equations
32
33 if err >0
34 %display the final result
35 fprintf(‘===================================\n’)
46 AntoineT = [6.95464, 1344.8, 219.48];
47 pEB = 250;
48 pT = 343;
49
50 %unpack z for convenience
61 AntoineEB = [6.95719, 1424.255, 213.21];
62 AntoineT = [6.95464, 1344.8, 219.48];
63 pEB = 250;
64 pT = 343;
65
76
77
78 function out = getPsat(antoine,T)
79 %handy function to get Psat
80 A = antoine(1);
90 C = antoine(3);
91 Psat = getPsat(antoine,T);
92 chain = B*log(10)/(T+C)ˆ2; %part from the chain rule
93 out = Psat*chain;
Note that I broke up the calculations of the saturation pressures and their deriva-
tives into separate functions. This is just for convenience, since it makes the program
easier to read. It is certainly not a requirement.
The program exhibits quadratic convergence to the answer:
Starting Newton-Raphson
===================================
k norm(R)
===================================
(3.27) The files for this problem are in the folder s13c5p123 matlab
(a) For a counter-current heat exchanger, using the notation from transport, the
equations we need are
We can use this equation to solve for the cooling fluid outlet temperature. To
compute the area, we use the design equation
The numerical calculations are in the Matlab file.
1function s13c5p1
2clc
3close all
4
8%T1 = inlet of center fluid
9%T2 = outlet of center fluid
10 %T1p = outlet of cooling fluid
11 %T2p = inlet of cooling fluid
12
23
24 %set outlet of center fluid
25 T2 = 50;
26
27 %first law balance to get outlet temperature
37 dTlm = (d2-d1)/log(d2/d1);
38 fprintf(‘The log-mean temperature difference is %5.3f …
C\n’,dTlm)
39 A = Q/U/dTlm;
40 fprintf(‘The area is %6.4f mˆ2\n’,A)
The output to the screen is
1The heat transfer is 345.0 kW
(b) To get a single non-linear equation, we need to eliminate T′
1using the equations
for Q,
Nonlinear equations 155
ln β+φT2−T′
T2−T′
2!
For Newton’s method, we also need the derivative
The quadratic convergence is seen in the last two steps in the m-file output, where
the error goes from 10−5to 10−14.
The Matlab file is:
1function s13c5p2
2clc
7%wp, Cpp = cooling fluid
8%T1 = inlet of center fluid
9%T2 = outlet of center fluid
20 %inlet parameters
21 T1 = 100;%deg C
22 T2p = 15;% deg C
23
24 %set outlet of center fluid
33 %compute the area
34 U = 1000; %W/m2 K
35 d2 = T2 – T2p;
36 d1 = T1 – T1p;
37 dTlm = (d2-d1)/log(d2/d1);
47 a=w*Cp/(U*A);
48 b = T1*(1-p);
49 %compute the temperature with Newton’s method
50 getT(T2,T1,T2p,p,a,b,w);
51
%6.4e\n’,count,T2,err)
59 %%%%%%%%%%%%%%%%%%
%6.4e\n’,count,T2,err);
63 %%%%%%%%%%%%%%%%%%
64 while abs(err) >1e-8
%6.4e\n’,count,T2,err)
%6.4e\n’,count,T2,err);
74 %%%%%%%%%%%%%%%%%%
Nonlinear equations 157
75 if count >20
76 fprintf(‘Did not converge! \n’)
87 out = a*(T1-T2)-num/denom;
88
89 function out = getdf(T2,T1,T2p,p,a,b)
90 c = log((b + p*T2 – T2p)/(T2-T2p))*(p – 1);
91 d = b-T2*(1-p);
6
7Starting Newton Method for w = 3.000e+00
11 k = 3 T2 = 49.99999361 err = 1.1259e-05
12 k = 4 T2 = 50.00000000 err = 9.9476e-14
(c) The Matlab program is:
1function s13c5p3
8%wp, Cpp = cooling fluid
9%T1 = inlet of center fluid
10 %T2 = outlet of center fluid
11 %T1p = outlet of cooling fluid
12 %T2p = inlet of cooling fluid
22 T1 = 100;%deg C
23 T2p = 15;% deg C
24
25 %set outlet of center fluid
26 T2 = 50;
36 d2 = T2 – T2p;
37 d1 = T1 – T1p;
38 dTlm = (d2-d1)/log(d2/d1);
39 fprintf(‘The log-mean temperature difference is %5.3f …
C\n’,dTlm)
49 %at very low flow rate, a good guess is a very slight …
change in the
50 %temperature of the inlet. so we should set the first guess …
to be
51 T2 = T1-0.1;