Assignment 3 - 2010 - Solution

From Process Model Formulation and Solution

Jump to: navigation, search

Question 1 [2]

The Secant method is the most widely used algorithm for solving a nonlinear equation, in chemical engineering, and other areas. Write a MATLAB or Python function that implements the Secant method, using a minimum number of function evaluations. (Python users: please do not just copy the built-in Secant method).

Your function should accept the following inputs:

  • a function handle to \(f(x)\) [see the software tutorial on the course website for a description of function handles]
  • a relaxation coefficient, \(0 < \alpha \leq 1.0\)
  • the stopping tolerance, \(\varepsilon_\text{tol}\)
  • the maximum number of iterations, maxiter (default = 200)
  • a single starting guess, \(x^{(0)}\)

Your code should automatically calculate \(x^{(-1)} = 1.001x^{(0)}\), and should use a minimum number of function evaluations. The stopping criterion should be as shown in the course notes:

\[\left|f(x^{(k)})\right| \leq \epsilon_{\rm tol} \qquad \text{and/or} \qquad \left|\frac{x^{(k)}-x^{(k-1)}}{x^{(k)}}\right| \leq \epsilon_{\rm tol}\]

Please print out a table at each iteration that shows these 5 outputs:

  • Iteration \(k\)
  • \(x^{(k)}\)
  • \(f(x^{(k)})\)
  • \(\displaystyle \left|\frac{x^{(k)}-x^{(k-1)}}{x^{(k)}}\right|\)
  • Relaxed \(\Delta x\)

Solution

See the following MATLAB and Python code:
MATLAB code Python code
%Uses the secant method to find f(x)=0.
function  x = secant_method(func, x0, alpha, tol, maxit)

%     INPUTS
%
%         f     : function f(x)
%         x0    : initial guess for x
%         alpha : relaxation coefficient: modifies Secant step size
%         tol   : convergence tolerance
%         maxit : maximum number of iteration, default=200

xprev = 1.001*x0;
x = x0;
f = @(x) func(x);

rel_step = 2.0 *tol;
k = 0;

fprintf('\n %10s %10s %12s %20s %20s\n','Iteration',' x',  ...
        'f(x)', 'Relative step', 'alpha * Delta x')

while (abs( f(x) ) > tol) && (rel_step) > tol && (k<maxit)
    
    rel_step = abs(x-xprev)/abs(x);
    
    % Full secant step
    dx = -f(x)/(f(x) - f(xprev))*(x - xprev);
    
    % Update `xprev` and `x` 
    xprev = x;
    x =  x + alpha*dx;
    
    fprintf('\n %5i %18f %12f %14f %19f\n',...
        k, xprev, f(xprev), rel_step, alpha*dx)
        
    k = k + 1;        
end
end
def secant_method(func, x0, alpha=1.0, tol=1E-9, maxit=200):
    """
    Uses the secant method to find f(x)=0.  
    
    INPUTS
    
        * f     : function f(x)
        * x0    : initial guess for x
        * alpha : relaxation coefficient: modifies Secant step size
        * tol   : convergence tolerance
        * maxit : maximum number of iteration, default=200        
    """
    # You can write more than one command per line in Python:
    x, xprev = x0, 1.001*x0
    f, fprev = func(x), func(xprev)
    
    rel_step = 2.0 *tol
    k = 0
    print('{0:12} {1:12} {2:12} {3:12} {4:12}'\
    .format('Iteration', 'x', 'f(x)', 'Rel step', 'alpha * Delta x'))
    
    while (abs(f) > tol) and (rel_step) > tol and (k<maxit):        
        rel_step = abs(x-xprev)/abs(x)
        
        # Full secant step
        dx = -f/(f - fprev)*(x - xprev)
        
        # Update `xprev` and `x` simultaneously
        xprev, x = x, x + alpha*dx
        
        # Update `fprev` and `f`:
        fprev, f = f, func(x)
        
        k += 1
        
        print('{0:10d} {1:12.5f} {2:12.5f} {3:12.5f} {4:12.5f}'\
        .format(k, xprev, fprev, rel_step, alpha*dx))
                
    return x

Question 2 [2]

In your reactor design course you will learn about adiabatic CSTR's exhibiting multiple steady states, depending on the reactor inlet temperature. In practice reactors are sometimes operated at an unstable steady point, due to process economics or improved conversion. The equation below was derived for such a system:

\[f(T) = (T-T_\text{feed}) - 150 \times \frac{1.34\times 10^9 e^{-\frac{7554}{T}}}{1+1.34\times 10^9 e^{-\frac{7554}{T}}}\]

There is more than one temperature at which \(f(T) = 0\); each \(f(T) = 0\) corresponds to a reactor steady state. Use the secant method, with relaxation factor \(\alpha = 0.9\), and \(T_\text{feed} = 300\) K to find one of these steady state temperatures, using starting values of \(T^{(0)}=320\) K and \(T^{(-1)}= 1.001 \times 320\) K.

  1. Show the first 3 iterations by hand. Explains what happens in the 3rd iteration, by means of a plot of \(f(x)\), and marking the iterations used so far on this plot.
  2. What could you change to avoid the problem in iteration 3?
  3. Compare your manual calculations in step 1 with the computer output from the previous question. Also report the temperature your computer code finds for \(f(T)=0\) at this starting guess. Select a reasonable tolerance value.
  4. Use the computer code from question 1, but with a starting value of \(T = 400\) K, and report the steady state operating temperature found.
  5. Comment on the solutions found from the two different starting points.

Solution

First we evaluate the function at the two initial guesses (given \(T_\text{feed}=300\)): \(f(T^{(0)}=320)=9.535\), and \(f(T^{(-1)})=1.001\times 320)= 9.623\). Now we apply the secant formula to find the first approximation \(T^{(1)}\):

\[T^{(1)}=T^{(0)}-\alpha\displaystyle\frac{f(T^{(0)})}{\displaystyle\frac{f(T^{(0)})-f(T^{(-1)})}{T^{(0)}-T^{(-1)}}}= 320.0 - 0.9 \displaystyle \frac{9.535}{\displaystyle \frac{9.535 - 9.623}{320 - 320.32}} = 320.0 - 0.9 \times 34.67 = 288.8\]

Evaluate \(f(T^{(1)})\), and repeat the above calculation to obtain \(T^{(1)}\):

\[\begin{split}T^{(2)}=T^{(1)}-\alpha\displaystyle\frac{f(T^{(1)})}{\displaystyle\frac{f(T^{(1)})-f(T^{(0)})}{T^{(1)}-T^{(0)}}}= 304.5\\\end{split}\]

Evaluate \(f(T^{(2)})\), and repeat the above calculation to obtain \(T^{(3)}\):

\[T^{(3)}=T^{(2)}-\alpha\displaystyle\frac{f(T^{(2)})}{\displaystyle\frac{f(T^{(2)})-f(T^{(1)})}{T^{(2)}-T^{(1)}}}=303.2\]
\[\]

The first three iterations (and the initial guess) are plotted below:

../_images/secant-iterations910.png

As you can see above, the first three iterations result in reasonable values, meaning that, as a sign of convergence, the difference between successive temperatures are getting less and less. Also, the value of the function itself is continually getting close to zero with the iterations. However, if you got, for example, negative absolute temperatures which is not physically acceptable, or the function value diverged from zero instead of getting closer to it, you could change your initial guesses or the \(\alpha\) value to avoid the problem.

The output from the code in question 1 is as follows:

Iteration    x            f(x)         Rel step     alpha * Delta x
         1    320.00000      9.53534      0.00100    -31.17468
         2    288.82532    -12.04971      0.10794     15.66275
         3    304.48807      1.18307      0.05144     -1.26028
         4    303.22778      0.24013      0.00416     -0.28886
         5    302.93892      0.02004      0.00095     -0.02368
         6    302.91525      0.00194      0.00008     -0.00228
         7    302.91296      0.00019      0.00001     -0.00023
         8    302.91274      0.00002      0.00000     -0.00002
         9    302.91271      0.00000      0.00000     -0.00000
        10    302.91271      0.00000      0.00000     -0.00000
302.912711645

The output from the code matches the hand calculations (otherwise something is wrong!). With the tolerance chosen as \(\varepsilon_\text{tol}=\sqrt{\text{eps}}\), the obtained temperature on convergence will be \(T^{(10)}=302.9\) K.

With a different initial guess as \(T^{(0)}=400.0\), and all other settings being the same, the obtained temperature will be: \(T^7=447.7\). You see that a different initial guess resulted in a far different solution. This is because the given function has multiple solutions. When an iterative method is used to solve the function, the obtained solution will depend on the choice of the initial guess. So, in order to find all solutions to a nonlinear equation, one way is to try different initial guesses. (Now try \(T^{(0)}=340.0\) and see what you get!)

The following computer code was used to solve this question:

MATLAB code Python code
% The function to be solved by the secant_method

f = inline(['(T - Tfeed) -150.*(1.34E9.*exp(-7554./T))' ...
             './(1+1.34E9.*exp(-7554./T))'], 'T', 'Tfeed');

Tfeed = 300;
maxiter = 200;
tol = sqrt(eps);
alpha = 0.9;
T0 = 320;
T = secant_method(@(T) f(T, Tfeed) ,T0,alpha,tol,maxiter)
import numpy as np
from matplotlib.pylab import *

# This is the code we wrote in question 1
# It was saved in a file called "secant.py"
from secant import secant_method

def f(T, Tfeed=300):
   return (T - Tfeed) - 150*(1.34E9*np.exp(-7554/T)) \
                         /(1+1.34E9*np.exp(-7554/T))

        
if __name__ == '__main__':
    
    # equivalent to sqrt(eps)
    tol = np.sqrt(np.finfo(np.double).eps)  
    root = secant_method(f, 320.0, 0.9, tol)
    print(root)
    
    # Create a plot
    fig = figure()
    T = np.arange(275, 350, 0.1)
    plot(T, f(T))
    grid('on')
    xlabel('Temperature [K]')
    ylabel('f(T)')
    
    # Plot iteration numbers 0, 1, 2, 3
    iter_T = np.array([320.0, 288.825, 304.488, 303.23])
    plot(iter_T, f(iter_T), 'ko', markersize=10 )
    for k, entry in enumerate([0, 1, 2, 3]):
        text(iter_T[k]-0.5, f(iter_T[k])-3.0, str(k))
    
    # Save the figure programatically 
    fig.savefig('../images/secant-iterations.png')

Question 3 [2]

In assignment 1 we derived the dynamic balance for a bioreactor system. The steady state nutrient concentration, \({\sf \overline{N}}\) and the steady state biomass concentration, \({\sf \overline{B}}\), can be shown to be:

\[\begin{split}\frac{dN}{dt} &= 0 = {\sf \overline{N}}_\text{in}\frac{F}{V} - {\sf \overline{N}}\frac{F}{V} - \left(\frac{1}{Y_{B}}\right)\left(\mu_{\text{max}}\frac{{\sf \overline{N}}}{K+{\sf \overline{N}}}\right){\sf \overline{B}}\\\frac{dB}{dt} &= 0 = - {\sf \overline{B}}\frac{F}{V} + \left(\mu_{\text{max}}\frac{{\sf \overline{N}}}{K+{\sf \overline{N}}}\right){\sf \overline{B}}\end{split}\]

Although we showed in assignment 1 that we can solve one equation and substitute it into the other, you should not do that for this question.

Rather, let \(\displaystyle \frac{dN}{dt} = f_1({\sf \overline{N}}, {\sf \overline{B}}) = 0\) and let \(\displaystyle \frac{dB}{dt} = f_2({\sf \overline{N}}, {\sf \overline{B}}) = 0\). Solve this coupled set of nonlinear equations for \({\bf x} = ({\sf \overline{N}}, {\sf \overline{B}})\) using Newton's method.

Use the following approach in your answer:

  1. Construct the Jacobian matrix, \({\bf J}\), symbolically in terms of the variables used in the equations.
  2. Simplify the Jacobian by substituting in these known constants:
    • \({\sf \overline{N}}_\text{in} = 100\) g/m3
    • \(F = 5000\) m3/day
    • \(V = 1600\) m3
    • \(Y_B = 0.80\) conversion efficiency
    • \(\mu_\text{max} = 5\) 1/day
    • \(K = 20\) g/m3
  3. Let your initial guess be \({\bf x}^{(0)} = (40 \text{g/m$^3$}, 20 \text{g/m$^3$})\). Calculate 2 iterations of the multivariate Newton-Raphson method, to show \({\bf x}^{(1)}\) and \({\bf x}^{(2)}\). Take a full Newton step at every iteration.
  4. Calculate and show the convergence tolerance, using the 2-norm for each iteration.
  5. Compare these iterations to the true solution presented by Elliot in assignment 1. Is Newton's method converging to this answer?

Feel free to code up the Jacobian matrix, and the functions \(f_1\) and \(f_2\) in MATLAB or Python and use your code to evaluate \({\bf J}({\sf \overline{N}}, {\sf \overline{B}})\), to avoid manual calculation errors.

Solution

  1. We start by deriving the Jacobian matrix for our two function system:

    \[\begin{split}J\left({\sf \overline{N}},{\sf \overline{B}}\right) &= \left[\begin{array}{ll} \frac{\partial f_1({\sf \overline{N}}, {\sf \overline{B}})}{\partial {\sf \overline{N}}} & \frac{\partial f_1({\sf \overline{N}}, {\sf \overline{B}})}{\partial {\sf \overline{B}}} \\ \frac{\partial f_2({\sf \overline{N}}, {\sf \overline{B}})}{\partial {\sf \overline{N}}} & \frac{\partial f_2({\sf \overline{N}}, {\sf \overline{B}})}{\partial {\sf \overline{B}}} \\\end{array}\right]\\J\left({\sf \overline{N}},{\sf \overline{B}}\right) &= \left[\begin{array}{ccc} - \frac{F}{V} - \left(\frac{1}{Y_{B}}\right)\left(\mu_{\text{max}}\frac{1}{K+{\sf \overline{N}}}\right){\sf \overline{B}} + \left(\frac{1}{Y_{B}}\right)\left(\mu_{\text{max}}\frac{{\sf \overline{N}}}{\left(K+{\sf \overline{N}}\right)^{2}}\right){\sf \overline{B}} & \quad , \quad & - \left(\frac{1}{Y_{B}}\right)\left(\mu_{\text{max}}\frac{{\sf \overline{N}}}{K+{\sf \overline{N}}}\right) \\ \left(\mu_{\text{max}}\frac{1}{K+{\sf \overline{N}}}\right){\sf \overline{B}} - \left(\mu_{\text{max}}\frac{{\sf \overline{N}}}{\left(K+{\sf \overline{N}}\right)^{2}}\right){\sf \overline{B}} & \quad , \quad & - \frac{F}{V} + \left(\mu_{\text{max}}\frac{{\sf \overline{N}}}{K+{\sf \overline{N}}}\right) \\\end{array}\right]\end{split}\]
  2. Subbing in the constant values provided for this question results in the following simplified Jacobian.

    \[\begin{split}J\left({\sf \overline{N}},{\sf \overline{B}}\right) &= \left[\begin{array}{ccc} - \frac{\left(5000 \frac{m^{3}}{day}\right)}{\left(1600 m^{3}\right)} - \left(\frac{1}{\left(0.80\right)}\right)\left(\left(5 \frac{1}{day}\right)\frac{1}{\left(20 \frac{gCOD}{m^{3}}\right)+{\sf \overline{N}}}\right){\sf \overline{B}} + \left(\frac{1}{\left(0.80\right)}\right)\left(\left(5 \frac{1}{day}\right)\frac{{\sf \overline{N}}}{\left(\left(20 \frac{gCOD}{m^{3}}\right)+{\sf \overline{N}}\right)^{2}}\right){\sf \overline{B}} & \quad , \quad & - \left(\frac{1}{\left(0.80\right)}\right)\left(\left(5 \frac{1}{day}\right)\frac{{\sf \overline{N}}}{\left(20 \frac{gCOD}{m^{3}}\right)+{\sf \overline{N}}}\right) \\ \left(\left(5 \frac{1}{day}\right)\frac{1}{\left(20 \frac{gCOD}{m^{3}}\right)+{\sf \overline{N}}}\right){\sf \overline{B}} - \left(\left(5 \frac{1}{day}\right)\frac{{\sf \overline{N}}}{\left(\left(20 \frac{gCOD}{m^{3}}\right)+{\sf \overline{N}}\right)^{2}}\right){\sf \overline{B}} & \quad , \quad & - \frac{\left(5000 \frac{m^{3}}{day}\right)}{\left(1600 m^{3}\right)} + \left(\left(5 \frac{1}{day}\right)\frac{{\sf \overline{N}}}{\left(20 \frac{gCOD}{m^{3}}\right)+{\sf \overline{N}}}\right) \\\end{array}\right]\end{split}\]
  3. In this case we are in the luckiest of situations, we know the exact solution so we can track both relative and absolute errors. The exact solution (taken from Assignment 1) is given by the following equations:

    \[\begin{split}{\sf \overline{N}} &= \left(\frac{\frac{F}{V_{R}}K}{\left(\mu_{\text{max}}-\frac{F}{V_{R}}\right)}\right)\\{\sf \overline{B}} &= \frac{\frac{F}{V_{R}}\left(N_{0} - {\sf \overline{N}}\right)}{\left(\frac{1}{Y_{B}}\right)\left(\mu_{\text{max}}\frac{{\sf \overline{N}}}{K+{\sf \overline{N}}}\right)}\end{split}\]

    Subbing in the given constants

    \[\begin{split}{\sf \overline{N}} &= \left(\frac{\left(\frac{\left(5000 \frac{m^{3}}{day}\right)}{\left(1600 m^{3}\right)}\right)\left(20 \frac{gCOD}{m^{3}}\right)}{\left(\left(5 \frac{1}{day}\right)-\left(\frac{\left(5000 \frac{m^{3}}{day}\right)}{\left(1600 m^{3}\right)}\right)\right)}\right)\\{\sf \overline{N}} &= 33.\overline{3} \frac{gCOD}{m^{3}}\\& \\{\sf \overline{B}} &= \frac{\left(\frac{\left(5000 \frac{m^{3}}{day}\right)}{\left(1600 m^{3}\right)}\right)\left(\left(100 \frac{gCOD}{m^{3}}\right) - \left(33.\overline{3} \frac{gCOD}{m^{3}}\right)\right)}{\left(\frac{1}{\left(0.80\right)}\right)\left(\left(5 \frac{1}{day}\right)\frac{\left(33.\overline{3} \frac{gCOD}{m^{3}}\right)}{\left(20 \frac{gCOD}{m^{3}}\right)+\left(33.\overline{3} \frac{gCOD}{m^{3}}\right)}\right)}\\{\sf \overline{B}} &= 53.\overline{3} \frac{gCOD}{m^{3}}\end{split}\]

    ITERATION 1

    Calculate Jacobian

    \[\begin{split}J^{(0)}\left(40 \frac{gCOD}{m^{3}},20 \frac{gCOD}{m^{3}}\right) &= \left[\begin{array}{ccc} -3.8194 \frac{1}{day} & \quad , \quad & -4.1667 \frac{1}{day}\\ 0.5556 \frac{1}{day} & \quad , \quad & 0.2083 \frac{1}{day}\\\end{array}\right]\end{split}\]

    Calculate function values

    \[\begin{split}f^{(0)}\left(40 \frac{gCOD}{m^{3}},20 \frac{gCOD}{m^{3}}\right) &= \left[\begin{array}{c} \frac{\left(5000 \frac{m^{3}}{day}\right)}{\left(1600 m^{3}\right)}\left(100 \frac{gCOD}{m^{3}}\right) - \frac{\left(5000 \frac{m^{3}}{day}\right)}{\left(1600 m^{3}\right)}\left(40 \frac{gCOD}{m^{3}}\right) - \left(\frac{1}{\left(0.80\right)}\right)\left(\left(5 \frac{1}{day}\right)\frac{\left(40 \frac{gCOD}{m^{3}}\right)}{\left(20 \frac{gCOD}{m^{3}}\right)+\left(40 \frac{gCOD}{m^{3}}\right)}\right)\left(20 \frac{gCOD}{m^{3}}\right) \\ - \frac{\left(5000 \frac{m^{3}}{day}\right)}{\left(1600 m^{3}\right)}\left(20 \frac{gCOD}{m^{3}}\right) + \left(\left(5 \frac{1}{day}\right)\frac{\left(40 \frac{gCOD}{m^{3}}\right)}{\left(20 \frac{gCOD}{m^{3}}\right)+\left(40 \frac{gCOD}{m^{3}}\right)}\right)\left(20 \frac{gCOD}{m^{3}}\right) \\\end{array}\right]\\&= \left[\begin{array}{c} 104.1667 \frac{gCOD/m^{3}}{day} \\ 4.1667 \frac{gCOD/m^{3}}{day} \\\end{array}\right]\end{split}\]

    Solve the system \(J\Delta x = -f\) (Use MATLAB or Python)

    \[\begin{split}\left[\begin{array}{ccc} -3.8194 \frac{1}{day} & \quad , \quad & -4.1667 \frac{1}{day}\\ 0.5556 \frac{1}{day} & \quad , \quad & 0.2083 \frac{1}{day}\\\end{array}\right] \Delta x &= -\left[\begin{array}{c} 104.1667 \frac{gCOD/m^{3}}{day} \\ 4.1667 \frac{gCOD/m^{3}}{day} \\\end{array}\right]\\\Delta x &= \left[\begin{array}{c} -25.7143 \frac{gCOD}{m^{3}} \\ 48.5714 \frac{gCOD}{m^{3}} \\\end{array}\right]\end{split}\]

    calculate new x

    \[\begin{split}x^{(1)} &= x^{(0)} + \Delta x \\x^{(1)} &= \left[\begin{array}{c} 40 \frac{gCOD}{m^{3}} \\ 20 \frac{gCOD}{m^{3}} \\\end{array}\right] + \left[\begin{array}{c} -25.7143 \frac{gCOD}{m^{3}} \\ 48.5714 \frac{gCOD}{m^{3}} \\\end{array}\right]\\x^{(1)} &= \left[\begin{array}{c} 14.2857 \frac{gCOD}{m^{3}} \\ 68.5714 \frac{gCOD}{m^{3}} \\\end{array}\right]\end{split}\]

    Calculate relative, absolute, and function error norms

    We will assume an \(\epsilon_{tol}\) of \(10^{-8}\) and for this solution we will use the 2-norm.

    Relative error

    \[\begin{split}\frac{\parallel x^{(1)} - x^{(0)} \parallel}{\parallel x^{(1)} \parallel} & < \epsilon_{tol} \; ? \\\frac{\parallel \Delta x \parallel}{\parallel x^{(1)} \parallel} & < \epsilon_{tol} \; ? \\\sqrt{\frac{\sum_{i=1}^{n} \left[ \Delta x_{i} \right]^{2}}{\sum_{i=1}^{n} \left[ x_{i}^{(1)} \right]^{2}}} & < \epsilon_{tol} \; ? \\\sqrt{\frac{(-25.7143)^{2} + (48.5714)^{2}}{(14.2857)^{2} + (68.5714)^{2}}} & < \epsilon_{tol} \; ? \\0.7846 & < \epsilon_{tol} \; ? \; \text{No, so keep going!}\end{split}\]

    Absolute error

    \[\begin{split}\frac{\parallel x^{(1)} - x^{(true)} \parallel}{\parallel x^{(true)} \parallel} & < \epsilon_{tol} \; ? \\\sqrt{\frac{\sum_{i=1}^{n} \left[ x_{i}^{(1)} - x_{i}^{(true)} \right]^{2}}{\sum_{i=1}^{n} \left[ x_{i}^{(true)} \right]^{2}}} & < \epsilon_{tol} \; ? \\\sqrt{\frac{(14.2857 - 33.\overline{3})^{2} + (68.5714 - 53.\overline{3})^{2}}{(33.\overline{3})^{2} + (53.\overline{3})^{2}}} & < \epsilon_{tol} \; ? \\0.3878 &< \epsilon_{tol} \; ? \; \text{No, so keep going!}\end{split}\]

    Function error

    \[\begin{split}\parallel f(x^{(1)}) \parallel & < \epsilon_{tol} \; ? \\\sqrt{\sum_{i=1}^{n} \left[ f(x^{(1)})_{i} \right]^{2}} & < \epsilon_{tol} \; ? \\\sqrt{(89.2857)^{2} + (-71.4286)^{2}} & < \epsilon_{tol} \; ? \\114.3415 &< \epsilon_{tol} \; ? \; \text{No, so keep going!}\end{split}\]

    ITERATION 2

    Calculate Jacobian

    \[\begin{split}J^{(1)}\left(14.2857 \frac{gCOD}{m^{3}},68.5714 \frac{gCOD}{m^{3}}\right) &= \left[\begin{array}{ccc} -10.4167 \frac{1}{day} & \quad , \quad & -2.6042 \frac{1}{day}\\ 5.8333 \frac{1}{day} & \quad , \quad & -1.0417 \frac{1}{day}\\\end{array}\right]\end{split}\]

    Calculate function values

    \[\begin{split}f^{(1)}\left(14.2857 \frac{gCOD}{m^{3}},68.5714 \frac{gCOD}{m^{3}}\right) &= \left[\begin{array}{c} 89.2857 \frac{gCOD/m^{3}}{day} \\ -71.4286 \frac{gCOD/m^{3}}{day} \\\end{array}\right]\end{split}\]

    Solve the system \(J\Delta x = -f\) (Use MATLAB or Python)

    \[\begin{split}\left[\begin{array}{ccc} -10.4167 \frac{1}{day} & \quad , \quad & -2.6042 \frac{1}{day}\\ 5.8333 \frac{1}{day} & \quad , \quad & -1.0417 \frac{1}{day}\\\end{array}\right] \Delta x &= -\left[\begin{array}{c} 89.2857 \frac{gCOD/m^{3}}{day} \\ -71.4286 \frac{gCOD/m^{3}}{day} \\\end{array}\right]\\\Delta x &= \left[\begin{array}{c} 10.7143 \frac{gCOD}{m^{3}} \\ -8.5714 \frac{gCOD}{m^{3}} \\\end{array}\right]\end{split}\]

    Calculate new x

    \[\begin{split}x^{(2)} &= x^{(1)} + \Delta x \\x^{(2)} &= \left[\begin{array}{c} 14.2857 \frac{gCOD}{m^{3}} \\ 68.5714 \frac{gCOD}{m^{3}} \\\end{array}\right] + \left[\begin{array}{c} 10.7143 \frac{gCOD}{m^{3}} \\ -8.5714 \frac{gCOD}{m^{3}} \\\end{array}\right]\\x^{(2)} &= \left[\begin{array}{c} 25.0000 \frac{gCOD}{m^{3}} \\ 60.0000 \frac{gCOD}{m^{3}} \\\end{array}\right]\end{split}\]

    Calculate relative, absolute, and function error norms

    We will assume an \(\epsilon_{tol}\) of \(10^{-8}\) and for this solution we will use the 2-norm.

    Relative error

    \[\begin{split}\frac{\parallel x^{(2)} - x^{(1)} \parallel}{\parallel x^{(2)} \parallel} & < \epsilon_{tol} \; ? \\\frac{\parallel \Delta x \parallel}{\parallel x^{(2)} \parallel} & < \epsilon_{tol} \; ? \\\sqrt{\frac{\sum_{i=1}^{n} \left[ \Delta x_{i} \right]^{2}}{\sum_{i=1}^{n} \left[ x_{i}^{(2)} \right]^{2}}} & < \epsilon_{tol} \; ? \\\sqrt{\frac{(10.7143)^{2} + (-8.5714)^{2}}{(25.0000)^{2} + (60.0000)^{2}}} & < \epsilon_{tol} \; ? \\0.2111 &< \epsilon_{tol} \; ? \; \text{No, so keep going!}\end{split}\]

    Absolute error

    \[\begin{split}\frac{\parallel x^{(2)} - x^{(true)} \parallel}{\parallel x^{(true)} \parallel} & < \epsilon_{tol} \; ? \\\sqrt{\frac{\sum_{i=1}^{n} \left[ x_{i}^{(2)} - x_{i}^{(true)} \right]^{2}}{\sum_{i=1}^{n} \left[ x_{i}^{(true)} \right]^{2}}} & < \epsilon_{tol} \; ? \\\sqrt{\frac{(25.0000 - 33.\overline{3})^{2} + (60.0000 - 53.\overline{3})^{2}}{(33.\overline{3})^{2} + (53.\overline{3})^{2}}} & < \epsilon_{tol} \; ? \\0.1697 &< \epsilon_{tol} \; ? \; \text{No, so keep going!}\end{split}\]

    Function error

    \[\begin{split}\parallel f(x^{(2)}) \parallel & < \epsilon_{tol} \; ? \\\sqrt{\sum_{i=1}^{n} \left[ f(x^{(2)})_{i} \right]^{2}} & < \epsilon_{tol} \; ? \\\sqrt{(26.0417)^{2} + (-20.8333)^{2}} & < \epsilon_{tol} \; ? \\33.3496 & < \epsilon_{tol} \; ? \; \text{No, so keep going!}\end{split}\]

    Note that all three errors appear to be decreasing from step to step so we are fairly confident the algorithm is converging and in fact this is true. If one were to continue the iterative process it would converge to within the specified relative tolerance within 7 iterations.

Personal tools
Namespaces
Variants
Actions
Navigation
Toolbox