Chapter 28, Problem 1P

Perform the first computation in Sec. 28.1, but for the casewhere $h=10$. Use the Heun (without iteration) and the fourth-orderRK method to obtain solutions.

To determine

To calculate: The solution for c if thedifferential equation for mass balance of single reactor is $V\frac{dc}{dt}=Q{c}_{in}-Qc$ by the Heun method and fourth-order RK method where $h=10$.

Answer to Problem 1P

Solution:

The solution for c by the Heun method where $h=10$ is,

Heun withoutiteration
t	c
0	10
10	25
20	34.375
30	40.23438
40	43.89648
50	46.1853

The solution for c by fourth-order RK method where $h=10$ is,

4th order RK
t	c
0	10
10	25.72917
20	35.27317
30	41.06419
40	44.57801
50	46.71009

Explanation of Solution

Given Information:

The differential equation for mass balance of single reactor is, $V\frac{dc}{dt}=Q{c}_{in}-Qc$.

The values,

$\begin{array}{l}V=100{\text{m}}^{\text{3}}\\ Q=5{\text{m}}^{\text{3}}/\text{min}\\ c_{in}=50\text{mg}/{\text{m}}^{\text{3}}\\ c_0=10\text{mg}/{\text{m}}^{\text{3}}\end{array}$

The analytical equation for mass balance of single reactor is,

$c=c_{in}\left(1-e^{-\left(Q/V\right)t}\right)+c_0{e}^{-\left(Q/V\right)t}$

Formula used:

The iteration formula for Heun’s method is,

$y_{i+1}=y_i+\frac{h}{2}\left(f\left(x_i,y_i\right)+f\left(x_i+h,y_i+hf\left(x_i,y_i\right)\right)\right)$

The fourth-order RK method for $\frac{dy}{dt}=f\left(t,y\right)$ is,

$\begin{array}{c}{y}_{n+1}=y_n+\frac{1}{6}\left(k_1+2k_2+2k_3+k_4\right)\\ t_{n+1}=t_n+h\end{array}$

Where,

$\begin{array}{c}{k}_1=hf\left(t_n,y_n\right)\\ k_2=hf\left(t_n+\frac{h}{2},y_n+\frac{k_1}{2}\right)\\ k_3=hf\left(t_n+\frac{h}{2},y_n+\frac{k_2}{2}\right)\\ k_4=hf\left(t_n+h,y_n+k_3\right)\end{array}$

Calculation:

Consider the analytical equation for mass balance of single reactor is,

$c=c_{in}\left(1-e^{-\left(Q/V\right)t}\right)+c_0{e}^{-\left(Q/V\right)t}$

Substitute the values $V=100{\text{ m}}^{\text{3}}\text{, }Q=5{\text{m}}^{\text{3}}/\text{min},c_{in}=50\text{mg}/{\text{m}}^{\text{3}}\text{ and }{c}_0=10\text{mg}/{\text{m}}^{\text{3}}$ in the above equation,

$\begin{array}{c}c=c_{in}\left(1-e^{-\left(5/100\right)t}\right)+c_0{e}^{-\left(5/100\right)t}\\ =50\left(1-e^{-0.05t}\right)+10e^{-0.05t}\end{array}$

Now, use VB code to determine c at different value of t using Heun’s method and RK4 method as below,

OptionExplicit

Subfind()

Dim t AsDouble, c AsDouble, h AsDouble,hhAsDouble

'Set the variables

t =0

c =10

h =10

'move to the cell b3

Range("b3").Select

ActiveCell. Value="Heun without iteration"

'Assign name to each columns

ActiveCell. Offset(1,0).Select

ActiveCell. Value="t"

ActiveCell. Offset(0,1).Select

ActiveCell. Value="c"

ActiveCell. Offset(0,1).Select

ActiveCell. Value="k_1"

ActiveCell. Offset(0,1).Select

ActiveCell. Value="c"

ActiveCell. Offset(0,1).Select

ActiveCell. Value="k_2"

ActiveCell. Offset(0,1).Select

ActiveCell. Value="fi"

'call Heun function to determine c at different values of t`

hh=Heun(t, c, h)

'Reset the values

t =0

c =10

h =10

'move to the cell b15

Range("b15").Select

'Assign name to each columns

ActiveCell. Value="4th order RK"

ActiveCell. Offset(1,0).Select

ActiveCell. Value="t"

ActiveCell. Offset(0,1).Select

ActiveCell. Value="c"

ActiveCell. Offset(0,1).Select

ActiveCell. Value="k_1"

ActiveCell. Offset(0,1).Select

ActiveCell. Value="cmid"

ActiveCell. Offset(0,1).Select

ActiveCell. Value="k_2"

ActiveCell. Offset(0,1).Select

ActiveCell. Value="cmid"

ActiveCell. Offset(0,1).Select

ActiveCell. Value="k_3"

ActiveCell. Offset(0,1).Select

ActiveCell. Value="cend"

ActiveCell. Offset(0,1).Select

ActiveCell. Value="k_4"

ActiveCell. Offset(0,1).Select

ActiveCell. Value="fi"

'call RK4 function to determine c at different values of t

hh= RK4(t, c, h)

EndSub

'Define the Heun function

FunctionHeun(t, c, h)

'Declare the variables

Dim dc1dt AsDouble,slpAsDouble, dc2dt AsDouble,ceAsDouble,cnewAsDouble

Dim j AsInteger

'Use loop to determine c at different value of t

For j =1To6

'move to cell b4

Range("b4").Select

'Display the value of t

ActiveCell. Offset(j,0).Select

ActiveCell. Value= t

'Display the value of c

ActiveCell. Offset(0,1).Select

ActiveCell. Value= c

'call drive function to determine derivative

dc1dt =drive(t, c)

ce= c +(dc1dt * h)

dc2dt =drive(t + h,ce)

slp=(dc1dt + dc2dt)/2

cnew= c +slp* h

t = t + h

'display the values in cell

ActiveCell. Offset(0,1).Select

ActiveCell. Value= dc1dt

ActiveCell. Offset(0,1).Select

ActiveCell. Value=ce

ActiveCell. Offset(0,1).Select

ActiveCell. Value= dc2dt

ActiveCell. Offset(0,1).Select

ActiveCell. Value=slp

c =cnew

Next

EndFunction

'define the drive functionto find the derivative

Functiondrive(t, c)

Dim Q AsDouble,cinAsDouble, v AsDouble, temp AsDouble

'set the variables

Q =5

cin=50

v =100

'use formula to find the derivative

temp =(Q *(cin- c))/ v

drive = temp

EndFunction

'define the RK4 function to find c

Function RK4(t, c, h)

Dim k_1 AsDouble, k_2 AsDouble, k_3 AsDouble, k_4 AsDouble

Dim cm AsDouble, cm1 AsDouble,ceAsDouble,slpAsDouble,cnewAsDouble

Dim j AsInteger

'determines k_1, k_2, k_3 and k_4 in Runge kutta to determine c

For j =1To6

'Move to cell b16

Range("b16").Select

ActiveCell. Offset(j,0).Select

ActiveCell. Value= t

ActiveCell. Offset(0,1).Select

ActiveCell. Value= c

'Call drive function to determine k_1

k_1 =drive(t, c)

cm = c +(k_1 *(h /2))

'Call drive function to determine k_2

k_2 =drive(t +(h /2), cm)

cm1 = c +(k_2 *(h /2))

'Call drive function to determine k_3

k_3 =drive(t +(h /2), cm1)

ce= c + k_3 * h

'Call drive function to determine k_4

k_4 =drive(t + h,ce)

slp=(k_1 +2*(k_2 + k_3)+ k_4)/6

cnew= c +(slp* h)

t = t + h

'Display values in cell

ActiveCell. Offset(0,1).Select

ActiveCell. Value= k_1

ActiveCell. Offset(0,1).Select

ActiveCell. Value= cm

ActiveCell. Offset(0,1).Select

ActiveCell. Value= k_2

ActiveCell. Offset(0,1).Select

ActiveCell. Value= cm1

ActiveCell. Offset(0,1).Select

ActiveCell. Value= k_3

ActiveCell. Offset(0,1).Select

ActiveCell. Value=ce

ActiveCell. Offset(0,1).Select

ActiveCell. Value= k_4

ActiveCell. Offset(0,1).Select

ActiveCell. Value=slp

c =cnew

Next

EndFunction

The following output gets displayed in the excel after the execution of the above code:

To draw the graph, use excel as below,

Step 1: Select cells from B5 to B10 and C5 to C10, then go to Insert tab and select the Line option from Charts subgroup.

Step 2: Select cells from B17 to B22 and C17 to C22, then go to Insert tab and select the Line option from Charts subgroup

Step 3: Merge the graphs.

The graph obtained is,

Hence, both the method gives the same results.