diff --git a/tutorial/julia/Tutorial 2 - Nonlinear Algebraic Systems/CSTR_ChloroBenzene.png b/tutorial/julia/Tutorial 2 - Nonlinear Algebraic Systems/CSTR_ChloroBenzene.png new file mode 100644 index 0000000..5767782 Binary files /dev/null and b/tutorial/julia/Tutorial 2 - Nonlinear Algebraic Systems/CSTR_ChloroBenzene.png differ diff --git a/tutorial/julia/Tutorial 2 - Nonlinear Algebraic Systems/J2_NLAsystems.ipynb b/tutorial/julia/Tutorial 2 - Nonlinear Algebraic Systems/J2_NLAsystems.ipynb index 4f8a150..4f55068 100644 --- a/tutorial/julia/Tutorial 2 - Nonlinear Algebraic Systems/J2_NLAsystems.ipynb +++ b/tutorial/julia/Tutorial 2 - Nonlinear Algebraic Systems/J2_NLAsystems.ipynb @@ -10,22 +10,21 @@ "\n", "# Learning Objectives \n", "\n", - "### *Setting up a linear system*\n", - "- Be able to distinguish a linear system from a nonlinear system.\n", - "- Learn to write a linear system in a matrix-vector format.\n", - "- Input a linear system into a program\n", - "- Write applicable engineering models in a linear format.\n", - "\n", - "### *Characterizing linear systems*\n", - "- Determine the bandwidth(s) of the linear system.*\n", - "- Determine if a linear system has a unique solution.\n", - "- Will solving a linear system provide a meaningful result (is it well-posed).\n", - "\n", - "### *Solving a linear system*\n", - "- Solve a linear system by Gauss-Elimination\n", - "- Solve a linear system by LU factorizations.\n", - "- State when numerical errors occur may occur in Gauss-Elimination and how to resolve them.\\*\n", - "- Be able to state some conditions for the convergence of Jacobi and Gauss-Siedel methods.\\*\n", + "By the end of this module, students will be able to\n", + "\n", + "### *Solve single nonlinear equations*\n", + "- Solve a single nonlinear equation by Picard’s method.\n", + "- Solve a single nonlinear equation by Newton’s method.\n", + "- Determine whether solving a nonlinear equation by fixed-point method will converge.\\*\n", + "- State the rate of convergence.\\*\n", + "\n", + "### *Solve systems of nonlinear equations*\n", + "- Formulate residual equations.\n", + "- Input a residual of a system of nonlinear equations in a program.\n", + "- Solve a system of nonlinear equations by Newton-Raphson.\n", + "\n", + "### *Solve parametric nonlinear systems by continuation methods*\n", + "- Solve a parametric nonlinear system using zeroth-order continuation for initial guesses.\\*\n", "\n", "\\*Not covered in tutorial.\n", "\n", @@ -33,461 +32,550 @@ ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "# Relevance \n", + "# Solving Single Nonlinear Algebra Equations \n", + "\n", + "In this section, we will explore Picard’s method and Newton’s method, which are two traditional methods to systematically solve a single nonlinear algebra equation $f(x)=0$.\n", + "\n", + "### Picard's Method\n", + "Picard’s method, which is also known as the method of successive substitution, involves simply adding $x$ to both sides of $f(x)=0$. The iterative scheme is\n", "\n", - "It would be only a slight overstatement to say that linear algebra underlies all modern numerical algorithms to one degree or another. Even software which specifically addresses nonlinear or complicated forms will often make use of linear algebra in a myriad of subroutines. In essence, many complicated problems can be reduced to repeated formulating and solving linear systems. \n", + "\\\\[x^{\\left(k+1\\right)} =f\\left(x^{\\left(k\\right)} \\right)+x^{\\left(k\\right)}.\\\\]\n", "\n", - "
" + "It is very easy to implement, but it may be very slow to converge.\n", + "\n", + "**Example 1**: Use Picard's method to find the root of function $f(x)=e^{-x}-x$ using initial guess $x^{(0)}=1$. \n", + "\n", + "We begin by defining a callback functions (f_ex1) to evaluate $f(x)$." ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": 1, "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "f_ex1 (generic function with 1 method)" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "# Chlorination of Benzene \n", - "\n", - "Monochlorobenzene is a high boiling point solvent used in many industrial applications ranging from the manufacture of [dyes](https://patentimages.storage.googleapis.com/76/55/12/31756682f755c5/US6013811.pdf), [plastic and rubber processing](https://patentimages.storage.googleapis.com/66/c8/e1/69918f557fbd35/US8563079.pdf), [paints](https://patentimages.storage.googleapis.com/28/17/a7/a40e46f75806cf/US9062206.pdf), and [waxes](https://patents.google.com/patent/US3999960A/en). Approximately, [$1 billion](https://www.prnewswire.com/news-releases/global-chlorobenzene-industry-300939916.html) worth of this solvent are used per year. It is [produced by combining chlorine with benzene](https://onlinelibrary.wiley.com/doi/abs/10.1002/14356007.o06_o03) in the presence of a Lewis acid catalyst such as Ferric Chloride FeCl3 or Anhydrous Aluminum Chloride AlCl3. When benzene and chlorine are mixed two key reactions occur: (1) the chlorination of benzene and (2) the secondary undesirable chlorination of the resulting monochlorobenzene. Product purity is ensured by performing a series of separation steps after the initial reaction and reagant usage is reduced by including a recycle loop. \n", - "\n", - "
\n", - "Objective: A large customer order has been placed for monochlorobenzene that needs to be fulfilled on a tight timeline. The small chemical company you work for has occasionally produced this solvent in the past and has an existing 6m3 continuously stirred tank reactor (CSTR) that will need to be repurposed to produce 25kmol/h monochlorobenzene. Moreover, there isn't sufficient development time available to attempt to improve the processing conditions for the reactor or separators. As a result, you'll need to use the documented process conditions for the reactor and separators which means supplying to the first separator an effluent consisting of 90% benzene, 7% monochlorobenzene, and 3% dichlorobenzene by mole. We need to determine the recycle flow rates and reaction rates required needed to achieve the desired rate at which monochlorobenzene need to be produced. \n", - "
\n", - "\n", - "Denote benzene as (A), monochlorobenzene as (B), and dichlorobenzene as (C). The mole fraction of species A in stream S1 is denote by y1,A and the molar flowrate F [kmol/h] for stream S1 is denoted F1.\n", - "\n", - "- The feed stream is composed entirely of A and Cl2.\n", - "- An excess amount of Cl2 is provided and HCl removal can be assumed to be trivial. So it's reasonable exclude it from our design considerations. Moreover, this means the reactions which take place in the reactor vessel can be simplified to A → B (at rate r1 [kmol/h]) and B → C (at rate r2 [kmol/h]).\n", - "\n", - "\n", - "\n", - "# Variables and Mass Balances \n", - "\n", - "For streams containing only a pure substance, we need only keep track of the total flow rate. Enumerating each quantity in the model is given by $F_1$, $F_2$, $F_3$, $F_4$, $F_5$, $F_6$, $F_7$, $y_{3,A}$, $y_{3,B}$, $y_{3,C}$, $y_{4,B}$, $y_{4,C}$, $r_1$, $r_2$, and $\\nu$. Five of these are specified (e.g. these are **parameters** with fixed values). This leaves ten degrees of freedom associated with a variable. We'll now need to write out the equations which govern the system (physical laws, equipment specifications, and performance specifications). We'll start with writing the mass balances over each unit operation. \n", - "\n", - "Mixer (Overall Mass Balance): \n", - "\n", - "\\begin{align}\n", - " F_1 + F_7 - F_2 &= 0\n", - "\\end{align}\n", - "\n", - "Reactor (Individual Component Mass Balances):\n", - "\n", - "\\begin{align}\n", - " F_2 - r_1\\nu - y_{3,A}F_3 &= 0 \\\\\n", - " (r_1 - r_2)\\nu + y_{3,B}F_3 &= 0 \\\\\n", - " r_2\\nu - y_{3,C}F_3 &= 0\n", - "\\end{align}\n", - "\n", - "Separator #1 (Individual Component Mass Balances):\n", - "\n", - "\\begin{align}\n", - " F_3 - F_4 - F_7 &= 0 \\\\\n", - " y_{3,B}F_3 - y_{4,B}F_4 &= 0 \\\\\n", - " y_{3,C}F_3 - y_{4,C}F_4 &= 0\n", - "\\end{align}\n", - "\n", - "Separator #2 (Individual Component Mass Balances):\n", - "\n", - "\\begin{align}\n", - " F_4 - F_5 - F_6 &= 0 \\\\\n", - " y_{4,B}F_4 - F_5 &= 0\n", - "\\end{align}\n", - "\n", - "In order to write these equations in a matrix-vector form: we'll first need to associate each variable (excluding **parameters**) in the model with a component in a vector of variables, $\\mathbf{x}$. \n", - "\n", - "\\begin{align}\n", - " \\mathbf{x} = (F_1, F_2, F_3, F_4, F_6, F_7, y_{4,B}, y_{4,C}, r_1, r_2)\n", - "\\end{align}\n", - "\n", - "So, $x_2 = F_2$, $x_3 = F_3$ and so on until we have $x_{10} = r_2$. \n", - "\n", - "
\n", - "Now, take a minute to inspect the model and categorize the system of equations. \n", - "
\n", - "\n", - "The 2nd and 3rd equations for Separator #1 as well as the 2nd equation in Separator #2 contain the expressions $y_{4,B}F_4$ and $y_{4,C}F_4$. These terms consist of two variables being multiplied by one another (a bilinear term) which is one of simplest commonly encountered nonlinear term. All other expressions consist of addition, subtraction, and the multiplication of a parameter by a variable. So this system of equations is nonlinear but would be linear if we ommitted the expressions: $y_{4,B}F_4$ and $y_{4,C}F_4$.\n", - "\n", - "
\n", - "It's often useful to reformulate a model (algebraicly rearrange and potentially introduce new variables) to arrive at a more easily solvable form. \n", - "
\n", - "\n", - "We'll start by introducing two new variables: the molar flowrates $F_{4,B}$ and $F_{4,C}$ defined as by\n", - "\n", - "\\begin{align}\n", - " F_{4,B} = y_{4,B}F_4 \\\\\n", - " F_{4,C} = y_{4,C}F_4\n", - "\\end{align}\n", - "\n", - "We can then write an additional equation for the molar flowrates of S4:\n", - "\n", - "Molar flowrate expressions for S4:\n", - "\n", - "\\begin{align}\n", - " F_{4,B} + F_{4,C} &= F_4\n", - "\\end{align}\n", - "\n", - "Next, we'll introduce the variables: $F_{3,A}$, $F_{3,B}$, and $F_{3,C}$ defined as:\n", - "\n", - "\\begin{align}\n", - " F_{3,A} = y_{3,A}F_3 \\\\\n", - " F_{3,B} = y_{3,B}F_3 \\\\\n", - " F_{3,C} = y_{3,C}F_3\n", - "\\end{align}\n", - "\n", - "We can then write an additional equation for the molar flowrates of S3:\n", - "\n", - "Molar flowrate expressions for S3:\n", - "\n", - "\\begin{align}\n", - "F_{3,A} + F_{3,B} + F_{3,C} = F_3\n", - "\\end{align}\n", - "\n", - "We can then rearrange the mass balances for the Reactor, Separator #1, and Separator #2.\n", - "\n", - "Reactor (Individual Component Mass Balances):\n", - "\n", - "\\begin{align}\n", - " F_2 - r_1\\nu - F_{3,A} &= 0 \\\\\n", - " (r_1 - r_2)\\nu + F_{3,B} &= 0 \\\\\n", - " r_2\\nu - F_{3,C} &= 0\n", - "\\end{align}\n", - "\n", - "Reformulated Separator #1 Mass Balances:\n", - "\n", - "\\begin{align}\n", - " F_3 - F_4 - F_7 &= 0 \\\\\n", - " F_{3,B} - F_{4,B} &= 0 \\\\\n", - " F_{3,C} - F_{4,C} &= 0\n", - "\\end{align}\n", - "\n", - "Reformulated Separator #2 Mass Balances:\n", - "\n", - "\\begin{align}\n", - " F_4 - F_5 - F_6 &= 0 \\\\\n", - " F_{4,B} - F_5 &= 0\n", - "\\end{align}\n", - "\n", - "\\begin{align}\n", - " \\mathbf{x'} = (F_1, F_2, F_{3,A}, F_{3,B}, F_{3,C}, F_3, F_{4,B}, F_{4,C}, F_4, F_6, F_7, r_1, r_2)\n", - "\\end{align}\n" + "f_ex1(x) = exp(-x) - x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "
\n", - "INTERACTIVE! Now we enter the above linear equation, Mx' = b in matrix-vector form.\n", - "
" + "We now define a function to perform Picard's method and run it for this problem." ] }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 2, "metadata": {}, "outputs": [ { - "data": { - "text/plain": [ - "-0.07" - ] - }, - "execution_count": 19, - "metadata": {}, - "output_type": "execute_result" + "name": "stdout", + "output_type": "stream", + "text": [ + "The solution is x = 0.5671477142601192\n" + ] } ], "source": [ - "# Define the Mx' = b linear system\n", - "\n", - "\n", - "# Define storage for linear system\n", - "M = zeros(13,13) # Makes an 2d array of size 13-by-13. This array has 13 rows and 13 columns.\n", - "b = zeros(13) # Makes a column vector of size 10. A vector with 13 rows and 1 column.\n", - "\n", - "# Input parameters\n", - "y3A = 0.90 # effluent benzene concentration of reactor\n", - "y3B = 0.07 # effluent monochlorobenzene concentration of reactor\n", - "y3C = 0.03 # effluent dichlorobenzene concentration of reactor\n", - "V = 6 # reactor volume\n", - "F_5 = 25.0 # flow rate of monochlorobenzene required\n", - "\n", - "### MASS BALANCE (MIXER) ###\n", - "\n", - "# fills in the first row\n", - "M[1,1] = 1.0 # Set the matrix M's entry in the first row and first column to 1\n", - "M[1,2] = -1.0 # Set the matrix M's entry in the first row and second column to -1\n", - "M[1,11] = 1.0 # Set the matrix M's entry in the first row and eighth column to 1\n", - "\n", - "# b[1] and all other entries of M in row #1, so no need to change these values\n", - "\n", - "M[2,2] = 1.0 # reactor - equation #1 \n", - "M[2,3] = -1.0\n", - "M[2,12] = -V\n", - "\n", - "M[3,4] = 1.0 # reactor - equation #2 \n", - "M[3,12] = -V\n", - "M[3,13] = V\n", - "\n", - "M[4,5] = -1.0 # reactor - equation #3 \n", - "M[4,13] = V\n", - "\n", - "M[10,8] = 1.0 # stream 4 balance \n", - "M[10,9] = -1.0\n", - "M[10,7] = 1.0\n", + "\"\"\"\n", + "nonlinear_picard\n", + "\n", + "Solves for `f(x) = 0` via Picard iteration where `f` is the residual, `x` is an initial guess, \n", + "`kmax` is the maximum number of iterations, `tol` is the absolute convergence tolerance. This returns\n", + "a tuple which holds solution `x` and a boolean value `flag` which is true if the algorithm successfully\n", + "converged and false otherwise.\n", + "\"\"\"\n", + "function nonlinear_picard(f,x,kmax,tol)\n", + " f_val = f(x)\n", + " k = 0 \n", + " while abs(f_val) > tol\n", + " x += f_val # Update x\n", + " k += 1 # Update iteration counter\n", + " f_val = f(x) # Compute function at next value\n", + " if (k > kmax)\n", + " println(\"Maximum iterations exceeded.\")\n", + " break\n", + " end\n", + " end\n", + " flag = (k > kmax || x == Inf || x == -Inf) ? false : true\n", + " return x, flag\n", + "end\n", "\n", - "M[11,3] = 1.0 # stream 3 balance #1 \n", - "M[11,4] = 1.0\n", - "M[11,5] = 1.0\n", - "M[11,6] = -1.0\n", + "x0 = 1.0 # initial guess\n", + "kmax = 20 # max number of iterations allowed\n", + "tol = 1e-5 # convergence tolerance\n", + "x,flag = nonlinear_picard(f_ex1,x0,kmax,tol)\n", + "println(\"The solution is x = $(x)\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Newton's Method\n", "\n", - "M[12,3] = 1.0 # stream 3 balance #2\n", - "M[12,6] = -y3A\n", + "Newton’s method is the most common method for solving a single nonlinear equation. The iterative scheme is\n", "\n", - "M[13,4] = 1.0 # stream 4 balance #3\n", - "M[13,6] = -y3B\n", + "\\\\[x^{(k+1)}=x^{(k)}- \\frac{f(x^{(k)})}{f'(x^{(k)})}.\\\\]\n", "\n", - "### SEPARATOR #1 & 2 - MASS BALANCES ###\n", + "As indicated by the scheme, the derivative $f'(x)$ is required for iterations.\n", "\n", - "# FILL IN THE REST BELOW HERE IN ROWS 5 to 9 of M!!!!!!!!" + "In order to solve **Example 1** using Newton's method, we should derive the derivative term $f'(x)=-e^{-x}-1$." ] }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "\"This Matrix (M) is \"" + "df_ex1 (generic function with 1 method)" ] }, + "execution_count": 4, "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/plain": [ - "13×13 Array{Float64,2}:\n", - " 1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0\n", - " 0.0 1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.0 0.0\n", - " 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -6.0 6.0\n", - " 0.0 0.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 6.0\n", - " 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 -1.0 0.0 -1.0 0.0 0.0\n", - " 0.0 0.0 0.0 1.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0\n", - " 0.0 0.0 0.0 0.0 1.0 0.0 0.0 -1.0 0.0 0.0 0.0 0.0 0.0\n", - " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 -1.0 0.0 0.0 0.0\n", - " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0\n", - " 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.0 -1.0 0.0 0.0 0.0 0.0\n", - " 0.0 0.0 1.0 1.0 1.0 -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", - " 0.0 0.0 1.0 0.0 0.0 -0.9 0.0 0.0 0.0 0.0 0.0 0.0 0.0\n", - " 0.0 0.0 0.0 1.0 0.0 -0.07 0.0 0.0 0.0 0.0 0.0 0.0 0.0" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/plain": [ - "\"The r.h.s vector (b) is \"" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/plain": [ - "13-element Array{Float64,1}:\n", - " 0.0\n", - " 0.0\n", - " 0.0\n", - " 0.0\n", - " 0.0\n", - " 0.0\n", - " 0.0\n", - " 25.0\n", - " 25.0\n", - " 0.0\n", - " 0.0\n", - " 0.0\n", - " 0.0" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/plain": [ - "13-element Array{Float64,1}:\n", - " 83.33333333333337 \n", - " 833.3333333333339 \n", - " 750.0000000000006 \n", - " 58.33333333333337 \n", - " 25.0 \n", - " 833.3333333333339 \n", - " 58.33333333333337 \n", - " 25.0 \n", - " 83.33333333333337 \n", - " 58.33333333333337 \n", - " 750.0000000000006 \n", - " 13.888888888888896\n", - " 4.166666666666667" - ] - }, - "metadata": {}, - "output_type": "display_data" + "output_type": "execute_result" } ], "source": [ - "# Run this to display the input matrix and vector\n", - "display(\"This Matrix (M) is \") \n", - "display(M)\n", - "display(\"The r.h.s vector (b) is \") \n", - "display(b)\n", - "x = M\\b\n", - "display(x)" + "df_ex1(x) = -exp(-x) - 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "# Existence and Uniqueness of the Solution \n", - "\n", - "Now that the linear system has been formulated, let's check to see whether the system has a solution. For a square matrix M, M system is invertible\n", - "\n", - "
\n", - " INTERACTIVE! Make use of the function det to check see if matrix M is invertible. Does the linear equation Mx = b have a unique solution? \n", - "
" + "We now define a function to perform Newton's method and run it for this problem." ] }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 8, "metadata": {}, "outputs": [ { - "data": { - "text/plain": [ - "-1.0799999999999992" - ] - }, - "execution_count": 21, - "metadata": {}, - "output_type": "execute_result" + "name": "stdout", + "output_type": "stream", + "text": [ + "The solution is x = 0.567143285989123\n" + ] } ], "source": [ - "using LinearAlgebra: det\n", + "\"\"\"\n", + "nonlinear_newton\n", + "\n", + "Solves for `f(x) = 0` via Newton's method where `f` is the residual, `df` is the derivative of the residual,\n", + "`x` is an initial guess, `kmax` is the maximum number of iterations, and `tol` is the absolute convergence tolerance.\n", + "This returns a tuple which holds solution `x` and a boolean value `flag` which is true if the algorithm successfully\n", + "converged and false otherwise.\n", + "\"\"\"\n", + "function nonlinear_newton(f,df,x,kmax,tol)\n", + " fval = f(x)\n", + " k = 0\n", + " while abs(fval) > tol\n", + " dfval = df(x)\n", + " x -= fval/dfval\n", + " k += 1\n", + " fval = f(x)\n", + " if k > kmax\n", + " println(\"Maximum iterations exceeded.\")\n", + " break\n", + " end\n", + " end\n", + " flag = (k > kmax || x == Inf || x == -Inf) ? false : true\n", + " return x, flag\n", + "end\n", "\n", - "# FILL IN THE REST BELOW HERE\n", - "det(M)" + "x0 = 1.0 # initial guess\n", + "kmax = 20 # max number of iterations allowed\n", + "tol = 1e-5 # convergence tolerance\n", + "x,flag = nonlinear_newton(f_ex1,df_ex1,x0,kmax,tol)\n", + "println(\"The solution is x = $(x)\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "# Well-posed problems and the Condition Number \n", - "\n", - "While a linear system, `Mx = b`, may have a unique solution, we also seek to understand whether using a particular algorithm to solve `Mx = b` will yeild an accurate solution. We often don't know `M` or `b` exactly. As a consequence, solutions which vary greatly when `M` or `b` are slightly perturbed may be suspect. This is generally assessed by evaluating the **condition number** of a linear system and is defined as: $\\text{cond}(M) \\equiv ||M|| \\times ||M^{-1}||$. A derivation of this is given in section 2.3.4 of Dorfman and Daoutidis, Numerical Methods with Chemical Engineering Applications which arrives at the following inequality:\n", - "\n", - "$||\\Delta x || \\leq \\text{cond}(M)\\bigg(\\frac{|| \\Delta b || || x ||}{|| b ||}\\bigg)$\n", - "\n", - "For a fixed condition number, $\\text{cond}(M)$, solution $x$, and $b$ if we slightly perturb $b$ by $\\Delta b$ then $x$ may change by at most an amount proportional to the condition number. Another way of interpreting this is by applying the following rule of thumb: the condition number $\\kappa$ means that the method looses $\\log_{10}{\\kappa}$ of accuracy relative to rounding error.\n", - "\n", - "If an application leads to an ill-posed problem (values of **M** and **b** may be known to very high precision), there are a multitude of ways this may be dealt with. This most common and often most robust approach is to apply a [**preconditioner**](http://www.mathcs.emory.edu/~benzi/Web_papers/survey.pdf). That is we find an appropriate matrix **Y** and multiply both sides of the original linear system to create an equivalent new linear system, `(YM)x = (Yb)`, which has a lower condition number. We can then solve the equivalent modified system. " + "Netwon's method exhibits quadratic convergence, that is a huge advantage, as the convergence accelerates rapidly near the solution." ] }, { + "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "
\n", - "INTERACTIVE! Make use of the function cond to compute the condition number using the $L^2$-norm. Is the matrix A well-conditioned? \n", - "
" + "# Solving Systems of Nonlinear Algebra Equations \n", + "\n", + "In this section, we introduce the numerical methods to solve systems of nonlinear algebra equations. The general form is given by:\n", + "\n", + "\\\\[R_1(x_1,x_2,\\ldots,x_n)=0\\\\\n", + "R_2(x_1,x_2,\\ldots,x_n)=0\\\\\n", + "\\ \\ \\ \\ \\ \\ \\ \\vdots\\\\\n", + "R_n(x_1,x_2,\\ldots,x_n)=0\\\\]\n", + "\n", + "We need to be careful to write the equations in a form such that they are equal to zero in order to apply the fixed-point methods. The compact notation of the system of nonlinear equations is $\\bf R(\\bf x)=\\bf0$, with $\\bf x$ being the vector of unknowns as $\\mathbf x = (x_1,x_2,\\ldots,x_n)$, and \n", + "\n", + "\\\\[\\mathbf{R}=\\left\\lbrack \\begin{array}{c}\n", + "R_1 \\left(x_1 ,x_2 ,\\dots ,x_n \\right)\\\\\n", + "R_2 \\left(x_1 ,x_2 ,\\dots ,x_n \\right)\\\\\n", + "\\vdots \\\\\n", + "R_n \\left(x_1 ,x_2 ,\\dots ,x_n \\right)\n", + "\\end{array}\\right\\rbrack\\\\]\n", + "\n", + "a vector-valued function which we will refer to as the **residual**.\n", + "\n", + "### Newton-Raphson\n", + "A common method for solving nonlinear systems of equations is the Newton–Raphson method. It relies on the same idea as Newton’s method but is now generalized to n dimensions. The iterative scheme is given by\n", + "\n", + "\\\\[\\mathbf{J}(\\mathbf{x}^{(k)})(\\mathbf{x}^{(k+1)}-\\mathbf{x}^{(k)})=-\\mathbf{R}(\\mathbf{x}^{(k)}),\\\\]\n", + "\n", + "where the $\\mathbf J$ is the Jacobian matrix defined as\n", + "\n", + "\\\\[\\mathbf{J}=\\left\\lbrack \\begin{array}{ccc}\n", + "\\frac{\\partial R_1 }{\\partial x_1 } & \\dots & \\frac{\\partial R_1 }{\\partial x_n }\\\\\n", + "\\vdots & \\ddots & \\vdots \\\\\n", + "\\frac{\\partial R_n }{\\partial x_1 } & \\dots & \\frac{\\partial R_n }{\\partial x_n }\n", + "\\end{array}\\right\\rbrack,\\\\]\n", + "\n", + "that is to be evaluated at the value $\\mathbf{x}^{(k)}$. In component form, the elements of the Jacobian are the partial derivatives $J_{ij}=\\frac{\\partial R_i}{\\partial x_j}$.\n", + "It is common to use a shorthand notation $\\mathbf{\\delta}^{(k+1)}=\\mathbf{x}^{(k+1)}-\\mathbf{x}^{(k)}$ for the difference between the previous values and the new values of $\\mathbf x$. The iterative scheme then becomes the linear system\n", + "\n", + "\\\\[\\mathbf{J}(\\mathbf{x}^{(k)})(\\mathbf{\\delta}^{(k+1)})=-\\mathbf{R}(\\mathbf{x}^{(k)}),\\\\]\n", + "\n", + "\\\\[\\mathbf x^{(k+1)}:=\\mathbf x^{(k)}+\\mathbf{\\delta}^{(k+1)}.\\\\]\n", + "\n", + "We have reduced solving the system of nonlinear equations into an iterative method where, at each step, we need to solve a system of linear equations. Since $\\mathbf{x}^{(k)}$ changes at each time step, we need to solve the linear system for different right-hand-side constant vectors and Jacobian matrices.\n", + "\n", + "The implementation of this method goes through the following steps.\n", + "\n", + "1. Create a vector with the residuals.\n", + "2. Compute the elements in the Jacobian. \n", + "3. Pick an initial guess $\\mathbf{x}^{(0)}$.\n", + "4.Iterate until either $\\| \\mathbf R \\|$ and/or $\\|\\mathbf \\delta\\|$ are small.\n", + "\n", + "**Example 2:** Use Netwon-Raphson to find the root of the system of equations\n", + "\n", + "\\\\[f_1(x_1,x_2)=e^{-x_1}-x_2, \\\\\n", + "f_2(x_1,x_2)=x_1+x_2^2-3x_2, \\\\]\n", + "\n", + "using the initial guess $x_1^{(0)}=0$ and $x_2^{(0)}=0$.\n", + "\n", + "**Solution:**\n", + "The residual vector for this problem is\n", + "\n", + "\\\\[\\mathbf{R}=\\left\\lbrack \\begin{array}{c}\n", + "e^{-x_1 } -x_2 \\\\\n", + "x_1 +x_2^2 -3x_2 \n", + "\\end{array}\\right\\rbrack\\\\]\n", + "\n", + "If we take the logical choice for the vector of unknowns to be $\\mathbf x =(x_1,x_2)$, then the Jacobian is\n", + "\n", + "\\\\[\\mathbf{J}=\\left\\lbrack \\begin{array}{cc}\n", + "-e^{-x_1 } & -1\\\\\n", + "1 & {2x}_2 -3\n", + "\\end{array}\\right\\rbrack\\\\]\n", + "\n", + "We then define corresponding callback functions R_ex2! and J_ex2! which compute the residual and it's jacobian, respectively: " ] }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "1494.7023893088221" + "J_ex2! (generic function with 1 method)" ] }, - "execution_count": 22, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "using LinearAlgebra: cond\n", + "function R_ex2!(out, x)\n", + " out[1] = exp(-x[1]) - x[2]\n", + " out[2] = x[1] + x[2]^2 -3*x[2]\n", + " return nothing\n", + "end\n", "\n", - "# FILL IN THE REST BELOW HERE" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Solving the linear equation with Gauss Elimination " + "function J_ex2!(out, x)\n", + " out[1,1] = -exp(-x[1])\n", + " out[1,2] = -1\n", + " out[2,1] = 1\n", + " out[2,2] = 2*x[2] - 3\n", + " return nothing\n", + "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now let's write a quick script to solve Mx = b. We'll do this in three steps. Defining a function to perform forward elimination, a function to back substitution, and lastly a main function to perform each of the prior functions in sequence." + "Then, we implement the iterative algorithm by:" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 11, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The solution is (x1,x2) = (0.982752934497163,0.37427931266177783)\n" + ] + } + ], "source": [ - "function forward_elimination!(M,b,n)\n", + "using LinearAlgebra: norm\n", + "\n", + "function nonlinear_newton_raphson(R!,J!,x,kmax,tol)\n", " \n", - " for k = 1:(n - 1)\n", - " for i = (k + 1):n\n", - " m = M[i,k]/M[k,k]\n", - " for j = (k + 1):n\n", - " M[i,j] = M[i,j] - m*M[k,j]\n", - " end\n", - " b[i] = b[i] - m*b[k]\n", + " # create intermediate storage\n", + " nx = length(x) \n", + " Rout = zeros(nx)\n", + " Jout = zeros(nx,nx)\n", + " \n", + " R!(Rout, x)\n", + " k = 0\n", + " while norm(Rout) > tol\n", + " J!(Jout, x) # compute Jacobian\n", + " del = -Jout\\Rout # calculate step\n", + " x += del # perform step\n", + " k += 1\n", + " R!(Rout, x) # evaluate residual\n", + " if k > kmax\n", + " println(\"Maximum iterations exceeded.\")\n", + " break\n", " end\n", " end\n", " \n", - " return nothing\n", - "end" + " flag = (k > kmax || x == Inf || x == -Inf) ? false : true\n", + " return x, flag\n", + "end\n", + "\n", + "x0 = [0.0; 0.0] # initial guess\n", + "kmax = 100 # max number of iterations allowed\n", + "tol = 1e-5 # convergence tolerance\n", + "x, flag = nonlinear_newton_raphson(R_ex2!,J_ex2!,x0,kmax,tol)\n", + "println(\"The solution is (x1,x2) = ($(x[1]),$(x[2]))\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "
\n", - "Comment: Note that while M, b are used in prior cells these variables don't interfere in the definition of forward_elimination!. The variables listed in the argument of forward_elimination!, that is (M,x,b,n), are evaluated in local scope when forward_elimination! is called.\n", - "
" + "# Case Study: Reactor-Separator-Recycle Process for Chlorination of Benzene \n", + "\n", + "In this section, we consider a reactor-separator-recycle process for chlorination of benzene as illustrated in the following figure:\n", + "\n", + "\n", + "\n", + "In the continuous stirred-tank reactor (CSTR), benzene $\\text C_6 \\text H_6$ (A) reacts with chlorine $\\text{Cl}_2$ to yield monochlorobenzene $\\text C_6 \\text H_5 \\text{Cl}$ (B) and dichlorobenzene $\\text C_6 \\text H_4 \\text{Cl}_2$ (C). The chlorine $\\text{Cl}_2$ is in excess, thus the reactions can be simplified as:\n", + "\n", + "\\\\[A\\longrightarrow B \\\\]\n", + "\\\\[B\\longrightarrow C \\\\]\n", + "\n", + "The reaction rates are given by:\n", + "\n", + "\\\\[r_1=k_1 \\frac{y_{3,A}}{y_{3,A} \\hat v_A + y_{3,B} \\hat v_B + y_{3,C} \\hat v_C},\\\\\n", + "r_2 =k_2 \\frac{y_{3,B}}{y_{3,A} \\hat v_A + y_{3,B} \\hat v_B + y_{3,C} \\hat v_C}.\\\\]\n", + "\n", + "The reactor volume is $v=6 \\text m^3$, the reaction rate constants are $k_1=0.4\\text h^{-1}$ and $k_2 = 0.055 \\text h^{-1}$, and the specific volumes of each species are $\\hat v_A=8.937\\times10^{-2} \\text m^3/\\text{kmol}$, $\\hat v_B=1.018\\times10^{-1} \\text m^3/\\text{kmol}$, and $\\hat v_C=1.130\\times10^{-1}\\text m^3/\\text{kmol}$. \n", + "\n", + "The purporse is to produce $50 \\text{kmol}$ per hour in the outlet stream of the final separator ($F_5=50 \\text{kmol/h}$) and we now want to solve for **the mole fractions of each species** ($y_{3,A}$, $y_{3,B}$, and $y_{3,C}$) in the stream $F_3$.\n", + "\n", + "To solve this problem, we analyze the degrees of freedom for this system:\n", + "- Unknowns: $F_1,F_2,F_3,F_4,F_5,F_6,F_7,y_{3,A},y_{3,B},y_{3,C},y_{4,B},y_{4,C},k_1,k_2,\\hat v_A, \\hat v_B, \\hat v_C, v$\n", + "- Process specifications: $F_5,k_1,k_2,\\hat v_A, \\hat v_B, \\hat v_C, v$\n", + "- Degrees of freedom: 18 unkowns - 7 specifications = 11 degrees of freedom\n", + "\n", + "Thus, 11 independent equations are required for solving the system. Based on the conservation of mass, we can set up the system of equations:\n", + "\n", + "- Mixer:\n", + "\n", + "\\\\[ F_1+F_7=F_2 \\\\]\n", + "\n", + "- Reactor:\n", + "\n", + "\\\\[F_2-r_1v=y_{3,A} F_3\\\\\n", + "0=(r_1-r_2)v+y_{3,B} F_3\\\\\n", + "0=r_2v - y_{3,C}F_3\\\\]\n", + "\n", + "- Separator 1: \n", + "\n", + "\\\\[F_3=F_4+F_7\\\\\n", + "y_{3,B}F_3=y_{4,B}F_4\\\\\n", + "y_{3,C}F_3=y_{4,C}F_4\\\\]\n", + "\n", + "- Separator 2:\n", + "\n", + "\\\\[F_4=F_5+F_6\\\\\n", + "y_{4,B}F_4=F_5\\\\]\n", + "\n", + "- Other Relationships\n", + "\n", + "\\\\[y_{3,A}+y_{3,B}+y_{3,C}=1\\\\\n", + "y_{4,B}+y_{4,C}=1\\\\]\n", + "\n", + "Since the compositions are unknown and the reaction rate laws are complex, the system is nonlinear.\n", + "\n", + "### Solutions:\n", + "To solve a nonlinear system of equations, we can use the Newton-Raphson method. The 11 independent equations should be written in residual form by moving every term to one side:\n", + "\n", + "We define the variable vector for this problem as $\\mathbf x=(F_1,F_2,y_{3,A},y_{3,B},y_{3,C},F_3,y_{4,B},y_{4,C},F_4,F_6,F_7)$, then the residual vector is constructed as:\n", + "\n", + "\\\\[\\mathbf{R}\\left(\\mathbf{x}\\right)=\\left\\lbrack \\begin{array}{l}\n", + "x_1 -x_2 +x_{11} \\\\\n", + "x_2 -\\frac{k_1 x_3 v}{x_3 \\hat{\\;v_A } +x_4 \\hat{\\;v_B } +x_5 \\hat{\\;v_C } }-x_3 x_6 \\\\\n", + "\\frac{\\left({k_2 x_4 -k}_1 x_3 \\right)v}{x_3 \\hat{\\;v_A } +x_4 \\hat{\\;v_B } +x_5 \\hat{\\;v_C } }+x_4 x_6 \\\\\n", + "\\frac{k_2 x_4 v}{x_3 \\hat{\\;v_A } +x_4 \\hat{\\;v_B } +x_5 \\hat{\\;v_C } }+x_5 x_6 \\\\\n", + "x_6 -x_9 -x_{11} \\\\\n", + "x_4 x_6 -x_7 x_9 \\\\\n", + "x_5 x_6 -x_8 x_9 \\\\\n", + "x_9 -F_5 -x_{10} \\\\\n", + "x_7 x_9 -F_5 \\\\\n", + "x_3 +x_4 +x_5 -1\\\\\n", + "x_7 +x_8 -1\n", + "\\end{array}\\right\\rbrack =\\mathbf{0}\\\\]\n", + "\n", + "Then, we form the Jacobian as:\n", + "\n", + "\\\\[\\mathbf{J}\\left(\\mathbf{x}\\right)=\\left\\lbrack \\begin{array}{ccc}\n", + "\\frac{\\partial R_1 }{\\partial x_1 } & \\dots & \\frac{\\partial R_1 }{\\partial x_{11} }\\\\\n", + "\\vdots & \\ddots & \\vdots \\\\\n", + "\\frac{\\partial R_{11} }{\\partial x_1 } & \\dots & \\frac{\\partial R_{11} }{\\partial x_{11} }\n", + "\\end{array}\\right\\rbrack\\\\]\n", + "\n", + "Most entries in the Jacobian are zero, the non-zero entries are indicated in the code below. " ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": 18, "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "J_case! (generic function with 1 method)" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "
\n", - "INTERACTIVE! The overall gauss elimination function is defined below. You just need to finish coding the back_substitution! function.\n", - "
" + "function R_case!(out,x)\n", + " \n", + " # process specifications and constants\n", + " V = 6 # m^3\n", + " F5 = 25 # kmol/h\n", + " Va = 8.937e-2 # m^3/kmol\n", + " Vb = 1.018e-1 # m^3/kmol\n", + " Vc = 1.13e-1 # m^3/kmol\n", + " k1 = 0.4 # h^-1\n", + " k2 = 0.055 # h^-1\n", + "\n", + " out[1] = x[1] - x[2] + x[11] # (eqn 1) x1 - x2 + x11 = 0\n", + "\n", + " # (eqn 2) x2 - k1*x3*v/(x3*Va + x4*Vb + x5*Vc) - x3*x6 = 0\n", + " out[2] = x[2] - k1*x[3]*V/(x[3]*Va + x[4]*Vb + x[5]*Vc) - x[3]*x[6]\n", + "\n", + " # (eqn 3) (k2*x4 - k1*x3)*v/(x3*Va + x4*Vb + x5*Vc) + x4*x6 = 0\n", + " out[3] = (k2*x[4] - k1*x[3])*V/(x[3]*Va + x[4]*Vb + x[5]*Vc) + x[4]*x[6]\n", + "\n", + " # (eqn 4) k2*x4*v/(x3*Va + x4*Vb + x5*Vc) - x5*x6 = 0\n", + " out[4] = k2*x[4]*V/(x[3]*Va + x[4]*Vb + x[5]*Vc) - x[5]*x[6]\n", + "\n", + " out[5] = x[6] - x[9] - x[11] # (eqn 5) x6 - x9 - x11 = 0\n", + " out[6] = x[4]*x[6] - x[7]*x[9] # (eqn 6) x4*x6 - x7*x9 = 0\n", + " out[7] = x[5]*x[6] - x[8]*x[9] # (eqn 7) x5*x6 - x8*x9 = 0\n", + " out[8] = x[9] - F5 - x[10] # (eqn 8) x9 - F5 - x10 = 0\n", + " out[9] = x[7]*x[9] - F5 # (eqn 9) x7*x9 - F5 = 0\n", + " out[10] = x[3] + x[4] + x[5] - 1 # (eqn 10) x3 + x4 + x5 - 1 = 0\n", + " out[11] = x[7] + x[8] - 1 # (eqn 11) x7 + x8 - 1 = 0\n", + " \n", + " return nothing\n", + "end\n", + "\n", + "function J_case!(out,x)\n", + " \n", + " # process specifications and constants\n", + " V = 6 # m^3\n", + " F5 = 25 # kmol/h\n", + " Va = 8.937e-2 # m^3/kmol\n", + " Vb = 1.018e-1 # m^3/kmol\n", + " Vc = 1.13e-1 # m^3/kmol\n", + " k1 = 0.4 # h^-1\n", + " k2 = 0.055 # h^-1\n", + "\n", + " # Set all values in out to zero\n", + " fill!(out, 0.0)\n", + " \n", + " out[1,1] = 1; out[1,2] = -1; out[1,11] = 1 # dR1/dx1, dR1/dx2, and dR1/dx11\n", + "\n", + " # dR2/dx2, dR2/dx3, dR2/dx4, dR2/dx5, and dR2/dx\n", + " out[2,2] = 1\n", + " out[2,3] = -(k1*V*(x[3]*Va + x[4]*Vb + x[5]*Vc) - k1*x[3]*V*Va)/(x[3]*Va + x[4]*Vb + x[5]*Vc)^2 - x[6]\n", + " out[2,4] = k1*x[3]*V*Vb/(x[3]*Va + x[4]*Vb + x[5]*Vc)^2\n", + " out[2,5] = k1*x[3]*V*Vc/(x[3]*Va + x[4]*Vb + x[5]*Vc)^2\n", + " out[2,6] = -x[3]\n", + "\n", + " \n", + " # dR3/dx3, dR3/dx4, dR3/dx5, and dR3/dx6\n", + " out[3,3] = (-k1*V*(x[3]*Va + x[4]*Vb + x[5]*Vc) - (k2*x[4] - k1*x[3])*V*Va)/(x[3]*Va + x[4]*Vb + x[5]*Vc)^2\n", + " out[3,4] = (k2*V*(x[3]*Va + x[4]*Vb + x[5]*Vc) - (k2*x[4] - k1*x[3])*V*Vb)/(x[3]*Va + x[4]*Vb + x[5]*Vc)^2 + x[6]\n", + " out[3,5] = -(k2*x[4] - k1*x[3])*V*Vc/(x[3]*Va + x[4]*Vb + x[5]*Vc)^2\n", + " out[3,6] = x[4]\n", + " \n", + " # dR4/dx3, dR4/dx4, dR4/dx5, and dR4/dx6\n", + " out[4,3] = -k2*x[4]*V*Va/(x[3]*Va + x[4]*Vb + x[5]*Vc)^2\n", + " out[4,4] = (k2*V*(x[3]*Va + x[4]*Vb + x[5]*Vc) - k2*x[4]*V*Vb)/(x[3]*Va + x[4]*Vb + x[5]*Vc)^2\n", + " out[4,5] = -k2*x[4]*V*Vc/(x[3]*Va + x[4]*Vb + x[5]*Vc)^2 - x[6]\n", + " out[4,6] = -x[5]\n", + "\n", + " # dR5/dx6, dR5/dx9, and dR5/dx11\n", + " out[5,6] = 1; out[5,9] = -1; out[5,11] = -1\n", + "\n", + " # dR6/dx4, dR6/dx6, dR6/dx7, and dR6/dx9\n", + " out[6,4] = x[6]; out[6,6] = x[4]; out[6,7] = -x[9]; out[6,9] = -x[7]\n", + "\n", + " # dR7/dx5, dR7/dx6, dR7/dx8, and dR7/dx9\n", + " out[7,5] = x[6]; out[7,6] = x[5]; out[7,8] = -x[9]; out[7,9] = -x[8]\n", + "\n", + " out[8,9] = 1; out[8,10] = -1 # dR8/dx9, dR8/dx10\n", + " out[9,7] = x[9]; out[9,9] = x[7] # dR9/dx7, dR9/dx9\n", + " out[10,3] = 1; out[10,4] = 1; out[10,5] = 1 # dR10/dx3, dR10/dx4, dR10/dx5\n", + " out[11,7] = 1; out[11,8] = 1 # dR11/dx7 and dR11/dx8\n", + " \n", + " return nothing\n", + "end" ] }, { @@ -496,124 +584,57 @@ "metadata": {}, "outputs": [], "source": [ - "function back_substitution!(M,x,b,n)\n", - " \n", - " # FILL IN THE REST BELOW HERE\n", - " \n", - " return nothing\n", - "end" + "Then, we implement the Newton-Raphson algorithm by:" ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 19, "metadata": {}, "outputs": [ { - "ename": "UndefVarError", - "evalue": "UndefVarError: M not defined", + "ename": "LoadError", + "evalue": "UndefVarError: R not defined", "output_type": "error", "traceback": [ - "UndefVarError: M not defined", + "UndefVarError: R not defined", "", "Stacktrace:", - " [1] top-level scope at In[2]:10" + " [1] R_case!(::Array{Float64,1}, ::Array{Float64,1}) at .\\In[18]:12", + " [2] nonlinear_newton_raphson(::typeof(R_case!), ::typeof(J_case!), ::Array{Float64,1}, ::Int64, ::Float64) at .\\In[11]:10", + " [3] top-level scope at In[19]:4", + " [4] include_string(::Function, ::Module, ::String, ::String) at .\\loading.jl:1091" ] } ], "source": [ - "function Gauss_Elimination!(M, b)\n", - " n = length(b)\n", - " x = zeros(n)\n", - " \n", - " forward_elimination!(M,b,n)\n", - " back_substitution!(M,x,b,n) \n", - " \n", - " return x\n", - "end\n", - "\n", - "# Runs a Gauss elimination routine using Mx = b has inputs\n", - "x = Gauss_Elimination!(M, b)\n", - "\n", - "display(\"The solution via Gauss elimination is \"); \n", - "display(x)" + "x0 = [35.7; 357; 0.9; 0.07; 0.03; 357; 0.7; 0.3; 35.7; 10.7; 321.3]\n", + "tol = 1e-4\n", + "kmax = 100\n", + "x,flag = nonlinear_newton_raphson(R_case!,J_case!,x0,kmax,tol)\n", + "println(\"The mole fraction of each species are\")\n", + "println(\"y_{3,A} = $(x[3])\") \n", + "println(\"y_{3,B} = $(x[4])\")\n", + "println(\"y_{3,C} = $(x[5])\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "# Solving with an LU factorization \n", - "\n", - "
\n", - "INTERACTIVE! Run the following snippet of code to solve the linear system using LU factorization to solve the linear system.\n", - "
" - ] - }, - { - "cell_type": "code", - "execution_count": 29, - "metadata": {}, - "outputs": [ - { - "data": { - "text/plain": [ - "\"The solution via LU factorization is \"" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/plain": [ - "13-element Array{Float64,1}:\n", - " 83.33333333333337 \n", - " 833.3333333333339 \n", - " 750.0000000000006 \n", - " 58.33333333333337 \n", - " 24.999999999999996\n", - " 833.3333333333339 \n", - " 58.33333333333337 \n", - " 25.0 \n", - " 83.33333333333337 \n", - " 58.33333333333337 \n", - " 750.0000000000006 \n", - " 13.888888888888896\n", - " 4.166666666666666" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "using LinearAlgebra: lu\n", - "\n", - "lu_fact = lu(M) # performs an LU factorization of M. Note that pivot = Val(false)\n", - " # means that the LU factorization is performed without pivoting\n", - " # (otherwise a permutation matrix describing the pivots is computed as well)\n", + "
\n", "\n", - "y = lu_fact.L\\b # solve Ly = b for y\n", - "x = lu_fact.U\\y # solve Ux = y for x\n", + "# Questions for reflection \n", "\n", - "display(\"The solution via LU factorization is \"); \n", - "display(x)" + "- XXX" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, - "source": [ - "
\n", - "\n", - "# Questions for reflection \n", - "\n", - "- What varieties of problems may result in a banded matrix?\n", - "- What does it mean that an iterative method converged? \n", - "- When is an absolute or relative convergence criteria preferable?\n", - "- When is solving by one method versus another preferable (e.g. Banded Gauss-Elimination versus Gauss-Siedel)?" - ] + "outputs": [], + "source": [] } ], "metadata": { diff --git a/tutorial/julia/Tutorial 3 - ODEs, Initial Value Problems/CSTR_pic.png b/tutorial/julia/Tutorial 3 - ODEs, Initial Value Problems/CSTR_pic.png new file mode 100644 index 0000000..c6e5fb8 Binary files /dev/null and b/tutorial/julia/Tutorial 3 - ODEs, Initial Value Problems/CSTR_pic.png differ diff --git a/tutorial/julia/Tutorial 4 - ODEs, Boundary Value Problems/J4_ODEBVPsystems.ipynb b/tutorial/julia/Tutorial 4 - ODEs, Boundary Value Problems/J4_ODEBVPsystems.ipynb index f1d85a7..1abae71 100644 --- a/tutorial/julia/Tutorial 4 - ODEs, Boundary Value Problems/J4_ODEBVPsystems.ipynb +++ b/tutorial/julia/Tutorial 4 - ODEs, Boundary Value Problems/J4_ODEBVPsystems.ipynb @@ -251,8 +251,15 @@ "# Questions for reflection \n", "\n", "- How can you tell if the spatial discretization chosen for this example was small enough?\n", - "- In general, a number of higher-order spatial discretizations may be used when discretizing a problem (see [slide 9 here](https://www.mathematik.uni-dortmund.de/~kuzmin/cfdintro/lecture4.pdf))." + "- In general, a number of higher-order spatial discretizations may be used when discretizing a problem (see [slide 9 here](https://www.mathematik.uni-dortmund.de/~kuzmin/cfdintro/lecture4.pdf)). What potential trade-offs may occur when using a higher-order accuracy scheme to discretize the problem solved in the=is tutorial? " ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/tutorial/matlab/Tutorial 2 - Nonlinear Algebraic Systems/M2_NLAsystems.mlx b/tutorial/matlab/Tutorial 2 - Nonlinear Algebraic Systems/M2_NLAsystems.mlx index b136892..75a2a7b 100644 Binary files a/tutorial/matlab/Tutorial 2 - Nonlinear Algebraic Systems/M2_NLAsystems.mlx and b/tutorial/matlab/Tutorial 2 - Nonlinear Algebraic Systems/M2_NLAsystems.mlx differ