978-1107135116 Chapter 4 Part 1

Nonlinear equations 199

325 %this is a vapor

326 P = vdw(V,T,a,b);

327 else

328 %this is on the horizontal part

339 ylabel(‘$P$ (MPa)’,‘FontSize’,14)

340 top = strcat(‘Isotherm for $T = $’,num2str(T),‘K’);

341 title(top,‘FontSize’,14)

342 out = h;

343

351 %vdW parameters: a,b

352 %units: P = MPa, V = L/mol, T = K

353

354 Pplot = zeros(500,1);

355 Vlplot = Pplot;

366 Vlplot(2) = Vl;

367 Vgplot(2) = Vg;

368 Tplot(2) = T;

369

370 fprintf(‘Go down from 395 K to 365K in -1 K increments’)

380

381 fprintf(‘Go up to critical point in +1 K increments’)

382 %now go up to 438 K from the value at 395 K. Reinitialize guesses

383 T = 395; Psat = Pplot(2); Vl = Vlplot(2); Vg = Vgplot(2);

384 for i = 33:75

395 for i = 76:125

396 T = T+0.1;

397 [Vl,Vg,Psat] = equilibrium(Psat,T,Vl,Vg,a,b);

398 Pplot(i) = Psat;

399 Vlplot(i) = Vl;

410 Vlplot(i) = Vl;

411 Vgplot(i) = Vg;

412 Tplot(i) = T;

413 end

414

425

426 %make the PV plot

427 g=figure;

428 plot(Vlplot,Pplot,‘ob’,Vgplot,Pplot,‘ob’)

429 xlabel(‘$\underline{V}$ (L/mol)’,‘FontSize’,14)

439 xlabel(‘$T$ (K)’,‘FontSize’,14)

440 ylabel(‘$P$ (MPa)’,‘FontSize’,14)

441 title(‘Phase envelope’,‘FontSize’,14)

442 text(370,1.8,‘Liquid’)

443 text(420,1,‘Vapor’)

(a) You want to rewrite the van der Waals equation of state in the form f(V) =0,

In Newton’s method, you need the derivative too,

A decent guess is 0.5, and the solution converged to a molar volume of 0.507813

L/mol.

The file output for this part is:

0.5 1 1.5 2 2.5 3

0

1

2

3

4

5

6

V(L/m ol)

P(MPa)

vdw Equat ion for T=470K

(b) Starting the volume very close to bwill usually converge to the root for the

liquid. Looking at the isotherm, you can see that the liquid root is around 0.4.

You can see that Newton’s method converged to the liquid’s root. The calculation

0.5 1 1.5 2 2.5 3

0

1

2

3

4

5

6

V(L/m ol)

P(MPa)

vdw Equat ion for T=425K

0 1 2 3 4 5 6 7 8 9 10

0.2

0.25

0.3

0.35

0.4

0.45

New t on method iterat ion, k

Value of molar volume, Vk(L/mol)

Evolution of Ne w ton Method f or V( 0 ) =0.23

10−14

10−12

10−10

10−8

10−6

10−4

10−2

100

102

Newton me tho d it e ration, k

Err or , |f(Vk|

Evolution of Ne w t on Method f or V( 0 ) =0.23

C o m p u t e d E rr or

S c a l in g R e s u l t

Nonlinear equations 203

(c) This calculation fails massively. If I print out the first 18 iterations, I can see

what happened:

Computing volume for Vinit = 0.600000000

k V err

1 0.704145230 -3.335338e-02

3 1.804569869 -5.655852e-01

5 0.803762387 -1.738042e-02

7 0.696989028 -3.650306e-02

9 1.122542959 -1.538455e-01

11 0.768445361 -1.765214e-02

13 0.832381352 -2.070888e-02

15 0.787732853 -1.685907e-02

The calculation fails at iteration 21 because there is an overﬂow error. The plot of

point. The roots all converge at the critical temperature, so it becomes extremely

diﬃcult to find the volume because the roots and the critical points of the function

are all very close. It is easy to hit the critical point of the function and get sent to

an unphysical, unrecoverable solution.

The file outputs for this part are:

204 Nonlinear equations

0 5 10 15 20 25

−14

−12

−10

−8

−6

−4

−2

0

2

x 1058

New t on method iterat ion, k

Value of molar volume, Vk(L/mol)

Evolution of Ne w ton Method f or V( 0 ) =0.6

−5 −4 −3 −2 −1 0 1 2 3 4 5

−20

−15

−10

−5

0

5

10

15

20

V(L/mol)

P(MPa)

vdw Eq uation f or T=425K

(d) The outputs for this part has a plot of the isotherm and a guess of 1.3 MPa

for the saturation pressure. At that pressure, a good guess for the volumes is 0.4

Vg−1

Vl

Vl−b

which has a very simple derivative,

f′(Psat)=Vl−Vg

Nonlinear equations 205

Note that you are working in somewhat strange units, but a MPa L/mol is actually

a kJ/mol. At each iteration, you want to first update Psat by Newton’s method,

k)

1.404886 MPa, Vl=0.352873 L/mol, and Vg=1.522065 L/mol. The plot of the

isotherm and the original cubic equation is provided.

0.5 1 1.5 2 2.5 3

0

0.5

1

1.5

2

2.5

3

V(L/mol)

P(MPa)

Guess of Psat for T=395K

2

3

4

5

6

P(MPa)

I s ot he r m f or T=395K

E O S v al u e s

I s ot h e r m

to 0.1K increments. To get really close to the critical point, I had to go to 0.01 K

206 Nonlinear equations

(4.1) The equations required for implicit Euler are

y1,i+1=y1,i−hky1,i+1y2

2,i+1

y2,i+1=y2,i−2hky1,i+1y2

2,i+1

y=0.9596. From the forward Euler, we have

yi+1=1+(0.01)(22−23)=0.96

0.1 0.5583 0.6 0.5149

1 -8.9167 -3 -17

(4.3) For one step of the predictor corrector with h=0.01 we get a predictor

yp

1=10 +0.01[cos(π)−102sin(π)] =9.99

(4.4) We first compute all of the kvalues:

k1=ln(2) =0.693

(4.5) Explicit Euler is

yn+1=yn+h f (yn)

and predictor-corrector is

In order for the values of yn+1to be the same, we need

Simplifying gives

The simplest way for this to be true is if f(yn)=0, i.e. if the derivative is zero at this

point and the function is ﬂat. It is also possible to have a function f(y) where the step

size is such that the current value and the explicit Euler value with a step size hgive

the same slope, such as a periodic function.

(4.6) The predictor-correction solution is

y1=yp

2=3+17/4

2=12 +17

8=29

8

The error in the solution with h=1/4 is |4−29/8|=3/8. Since the error in this

scheme is quadratic, we would expect that

Initial value problems 223

Solving for h, we get

(4.7) (a)

(b) Implicit Euler requires newtons method or something similar to solve for yi+1

First put it in the form of a residual

then diﬀerentiate with respect to yi+1

xis the independent variable so it is known. Newton’s method for the kth iteration

is:

i+1−yi−h(y2(k)

i+1exp(−2xi+1))

The converged result is the new value of yi+1

(c) Predictor-corrector starts with an explicit euler step

After is a second step using yp

i+1

the new value of yi+1is the average of the two values

yi+1=yp

i+1+yc

i+1

2

(4.8) (a)

224 Initial value problems

(c) The equation for Implicit Euler is

this equation must be solved for yi+1. One simple way to do this is Newtons

method. The value of xi+1is known so first we write the above equation in the

form of a residual.

i+1,xi+1)

The converged result of Newtons method is now the new value of yi+1

(4.9) (a) Explicit Euler requires 5 time steps to evaluate y(2.5) the first is:

(4.10) The script to used to generate the solution to this problem can be found in 92p6p34 matlab.

The Matlab program for this problem is

Initial value problems 225

1function f92p6p34

2clc

3close all

4%initialize time steps

15 yee(2,1)=6;

16 for i=2:npts+1 %loop for EE

17 rke=getk(yee(:,i-1));

18 yee(:,i)=yee(:,i-1)+h*rke;

19 tee(i)=tee(i-1)+h;

integration\n’)

30 disp(c) %displays the values for the initial condition and the …

time steps required

31

32 %initilization for RK4

43 rk3=getk(yrk(:,i-1)+(h/2)*rk2);

44 rk4=getk(yrk(:,i-1)+h*rk3);

45 yrk(:,i)= yrk(:,i-1)+(h/6)*(rk1+2*rk2+2*rk3+rk4);

46 trk(i)=trk(i-1)+h;

47 end

56 rk=[trk;yrk];

57 for i=1:3

58 str1=L(i,:);

59 str2=num2str(rk(i,:));

60 d(i,:)=[str1,str2];

70 out(1,1)=dy1;

71 out(2,1)=dy2;

(a) Explicit Euler requires 4 time steps to reach t=2. The values of the intermediate

steps are

(b) For RK4 only one time step is needed but 4 function evaluations are computed

the four K vectors (listed y1;y2) are

(4.11) (a) Explicit Euler requires 4 steps to reach y(1) the first is:

The next four are

Initial value problems 227

(b) RK4 takes only one time step but requires four function evaluations. They are:

The final value is y(1) ≈0.8203

(4.12) (a) What diﬀerential equation is solved by this program?

(e) What does the variable crepresent?

The derivative of the residual (or function) for Newton’s method. It also OK if

they just write the equation that we are describing:

(4.13) The local stability analysis is

∂f

(4.14) (a) Rewriting the equations in matrix form:

d

dt “y1

y2#=“0 1

−2−3#” y1

y2#

The constant matrix has eigenvalues with corresponding vectors:

Initial value problems 229

(d) The stability criteria for explicit euler is hλ≤2. The value of the eigenvalue

(4.15) Convert these into a system of equations by defining y1=yand y2=y′

The eigenvalue equation is λ2+4y2λ−1=0. Using the quadratic equation, the

eigenvalues are λ=−2y2±q4y2

2+1. The discriminant is always positive, so the

eigenvalues are real. For forward Euler to be stable, we require h≤2/|λmax|. The

maximum eigenvalue depends on the sign of the slope of the function:





−2y′−p(2y′)2+1 (y′>0)

(4.16) First we must compute the Jacobian in order to find the eigenvalues of problem about

Now we evaluate the Jacobian at t=0 such that (y1,y2)=(.5,5)

The eigenvalues of this matrix are λ=−3.618 and λ=−1.382. The stability criteria

230 Initial value problems

for explicit euler is hλ≤2. The value of the eigenvalue used is the absolute value of

the largest eigenvalue thus:

Based on the values of the eigenvalues returned for this matrix about t=0 I would

not expect this system to be stiﬀ.

(4.17) This problem is the forward Euler integration of the diﬀerential equation

For non-linear ODEs, the approximation for the eigenvalue at some point ysis

In this problem, the time step of 0.05 is definitely smaller than the max step size.

(4.18) You need 4 coupled equations. The explanation is not required, but the problem you

need to solve is:





(4.19) Let the subscript on ydenote the order of the derivative. For the first equation, we get

d



y0

y1





=



y1

y2





,



y0

y1





=



2

1





y3

y0

y3

0