\n",
+ "INTERACTIVE! Run the above snippet of code while varying the initial guess to see how the solution responds. Does this always yield the same value? If so, why? For extreme guesses are any errors encountered? Take a moment and try to reflect on condition(s) under which these may occur.\n",
+ "
"
+ ]
+ },
+ {
"cell_type": "markdown",
"metadata": {},
"source": [
@@ -353,6 +360,14 @@
"source": [
"using LinearAlgebra: norm\n",
"\n",
+ "\"\"\"\n",
+ "nonlinear_newton_raphson\n",
+ "\n",
+ "Solves for `R(x) = 0` via Newton's method where `R!` is the residual, `J!` is the Jacobian 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_raphson(R!,J!,x,kmax,tol)\n",
" \n",
" # create intermediate storage\n",
@@ -386,6 +401,7 @@
]
},
{
+ "attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
@@ -405,9 +421,9 @@
"\\\\[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",
+ "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",
+ "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",
@@ -447,13 +463,29 @@
"### 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",
+ "\\\\[\n",
+ "\\begin{array}{l}\n",
+ "F_1 -F_2 +F_7 =0 \\\\\n",
+ "F_2 - k_1 y_{3,A}v / (y_{3,A} {\\hat{v} }_A +y_{3,B} {\\hat{v} }_B +y_{3,C} {\\hat{v} }_C) - y_{3,A} F_3 =0\\\\\n",
+ "v(k_2 y_{3,B} -k_1 y_{3,A}) / (y_{3,A} {\\hat{v} }_A +y_{3,B} {\\hat{v} }_B +y_{3,C} {\\hat{v} }_C) + y_{3,B} F_3 =0\\\\\n",
+ "k_2 y_{3,B}v / (y_{3,A} {\\hat{v} }_A +y_{3,B} {\\hat{v} }_B +y_{3,C} {\\hat{v} }_C)-y_{3,C} F_3 =0 \\\\\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",
+ "F_4 -F_5 -F_6 =0\\\\\n",
+ "y_{4,B} F_4 -F_5 =0\\\\\n",
+ "y_{3,A} +y_{3,B} +y_{3,C} -1=0\\\\\n",
+ "y_{4,B} +y_{4,C} -1=0\n",
+ "\\end{array}\n",
+ "\\\\]\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_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",
@@ -471,12 +503,16 @@
"\\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. "
+ "Most entries in the Jacobian are zero, the non-zero entries are indicated in the code below. \n",
+ "\n",
+ "
\n",
+ "INTERACTIVE! Fill in the remaining entries to complete the residual and Jacobian below.\n",
+ "
"
]
},
{
"cell_type": "code",
- "execution_count": 18,
+ "execution_count": 20,
"metadata": {},
"outputs": [
{
@@ -485,7 +521,7 @@
"J_case! (generic function with 1 method)"
]
},
- "execution_count": 18,
+ "execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
@@ -513,9 +549,10 @@
" # (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",
+ " ### SPECIFY THE COMPONENTS 5 TO 7 OF THE RESIDUAL HERE ###\n",
+ " ### SPECIFY THE COMPONENTS 5 TO 7 OF THE RESIDUAL HERE ###\n",
+ " ### SPECIFY THE COMPONENTS 5 TO 7 OF THE RESIDUAL HERE ###\n",
+ " \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",
@@ -560,14 +597,9 @@
" 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",
+ " ### SPECIFY THE JACOBIAN OF COMPONENTS 5 TO 7 OF THE RESIDUAL HERE ###\n",
+ " ### SPECIFY THE JACOBIAN OF COMPONENTS 5 TO 7 OF THE RESIDUAL HERE ###\n",
+ " ### SPECIFY THE JACOBIAN OF COMPONENTS 5 TO 7 OF THE RESIDUAL HERE ###\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",
@@ -579,31 +611,25 @@
]
},
{
- "cell_type": "code",
- "execution_count": null,
+ "cell_type": "markdown",
"metadata": {},
- "outputs": [],
"source": [
"Then, we implement the Newton-Raphson algorithm by:"
]
},
{
"cell_type": "code",
- "execution_count": 19,
+ "execution_count": 21,
"metadata": {},
"outputs": [
{
- "ename": "LoadError",
- "evalue": "UndefVarError: R not defined",
- "output_type": "error",
- "traceback": [
- "UndefVarError: R not defined",
- "",
- "Stacktrace:",
- " [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"
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "The mole fraction of each species are\n",
+ "y_{3,A} = 0.9454833705152079\n",
+ "y_{3,B} = 0.05408780736715945\n",
+ "y_{3,C} = 0.00042882211763266596\n"
]
}
],
@@ -624,17 +650,10 @@
"source": [
" \n",
"\n",
- "# Questions for reflection \n",
+ "# Question(s) for reflection \n",
"\n",
- "- XXX"
+ "- The Newton-Raphson algorithm makes use of derivative information to compute the direction of steepest descent by solving a linear system of equations. One common issue encountered when applying this approach is potentially overshooting the soluton $\\mathbf{x} = \\mathbf{x}^*$. How might one solve this issue?"
]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {},
- "outputs": [],
- "source": []
}
],
"metadata": {
diff --git a/tutorial/julia/Tutorial 2 - Nonlinear Algebraic Systems/J2_NLAsystems_solved.ipynb b/tutorial/julia/Tutorial 2 - Nonlinear Algebraic Systems/J2_NLAsystems_solved.ipynb
new file mode 100644
index 0000000..f9523af
--- /dev/null
+++ b/tutorial/julia/Tutorial 2 - Nonlinear Algebraic Systems/J2_NLAsystems_solved.ipynb
@@ -0,0 +1,678 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Linear Algebraic Systems \n",
+ "\n",
+ " \n",
+ "\n",
+ "# Learning Objectives \n",
+ "\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",
+ " "
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# 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",
+ "\\\\[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": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "f_ex1 (generic function with 1 method)"
+ ]
+ },
+ "execution_count": 2,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "f_ex1(x) = exp(-x) - x"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We now define a function to perform Picard's method and run it for this problem."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "The solution is x = 0.5671477142601192\n"
+ ]
+ }
+ ],
+ "source": [
+ "\"\"\"\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",
+ "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",
+ "Newton’s method is the most common method for solving a single nonlinear equation. The iterative scheme is\n",
+ "\n",
+ "\\\\[x^{(k+1)}=x^{(k)}- \\frac{f(x^{(k)})}{f'(x^{(k)})}.\\\\]\n",
+ "\n",
+ "As indicated by the scheme, the derivative $f'(x)$ is required for iterations.\n",
+ "\n",
+ "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": 4,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "df_ex1 (generic function with 1 method)"
+ ]
+ },
+ "execution_count": 4,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "df_ex1(x) = -exp(-x) - 1"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We now define a function to perform Newton's method and run it for this problem."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "The solution is x = 0.567143285989123\n"
+ ]
+ }
+ ],
+ "source": [
+ "\"\"\"\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",
+ "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": [
+ "Netwon's method exhibits quadratic convergence, that is a huge advantage, as the convergence accelerates rapidly near the solution."
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "
\n",
+ "INTERACTIVE! Run the above snippet of code while varying the initial guess to see how the solution responds. Does this always yield the same value? If so, why? For extreme guesses are any errors encountered? Take a moment and try to reflect on condition(s) under which these may occur.\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# 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": 9,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "J_ex2! (generic function with 1 method)"
+ ]
+ },
+ "execution_count": 9,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "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",
+ "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": [
+ "Then, we implement the iterative algorithm by:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "The solution is (x1,x2) = (0.982752934497163,0.37427931266177783)\n"
+ ]
+ }
+ ],
+ "source": [
+ "using LinearAlgebra: norm\n",
+ "\n",
+ "\"\"\"\n",
+ "nonlinear_newton_raphson\n",
+ "\n",
+ "Solves for `R(x) = 0` via Newton's method where `R!` is the residual, `J!` is the Jacobian 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_raphson(R!,J!,x,kmax,tol)\n",
+ " \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",
+ " 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]))\")"
+ ]
+ },
+ {
+ "attachments": {},
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# 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",
+ "\\\\[\n",
+ "\\begin{array}{l}\n",
+ "F_1 -F_2 +F_7 =0 \\\\\n",
+ "F_2 - k_1 y_{3,A}v / (y_{3,A} {\\hat{v} }_A +y_{3,B} {\\hat{v} }_B +y_{3,C} {\\hat{v} }_C) - y_{3,A} F_3 =0\\\\\n",
+ "v(k_2 y_{3,B} -k_1 y_{3,A}) / (y_{3,A} {\\hat{v} }_A +y_{3,B} {\\hat{v} }_B +y_{3,C} {\\hat{v} }_C) + y_{3,B} F_3 =0\\\\\n",
+ "k_2 y_{3,B}v / (y_{3,A} {\\hat{v} }_A +y_{3,B} {\\hat{v} }_B +y_{3,C} {\\hat{v} }_C)-y_{3,C} F_3 =0 \\\\\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",
+ "F_4 -F_5 -F_6 =0\\\\\n",
+ "y_{4,B} F_4 -F_5 =0\\\\\n",
+ "y_{3,A} +y_{3,B} +y_{3,C} -1=0\\\\\n",
+ "y_{4,B} +y_{4,C} -1=0\n",
+ "\\end{array}\n",
+ "\\\\]\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. \n",
+ "\n",
+ "
\n",
+ "INTERACTIVE! Fill in the remaining entries to complete the residual and Jacobian below.\n",
+ "
"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "J_case! (generic function with 1 method)"
+ ]
+ },
+ "execution_count": 20,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "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"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Then, we implement the Newton-Raphson algorithm by:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "The mole fraction of each species are\n",
+ "y_{3,A} = 0.9454833705152079\n",
+ "y_{3,B} = 0.05408780736715945\n",
+ "y_{3,C} = 0.00042882211763266596\n"
+ ]
+ }
+ ],
+ "source": [
+ "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": [
+ " \n",
+ "\n",
+ "# Question(s) for reflection \n",
+ "\n",
+ "- The Newton-Raphson algorithm makes use of derivative information to compute the direction of steepest descent by solving a linear system of equations. One common issue encountered when applying this approach is potentially overshooting the soluton $\\mathbf{x} = \\mathbf{x}^*$. How might one solve this issue?"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Julia 1.5.1",
+ "language": "julia",
+ "name": "julia-1.5"
+ },
+ "language_info": {
+ "file_extension": ".jl",
+ "mimetype": "application/julia",
+ "name": "julia",
+ "version": "1.5.1"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/tutorial/julia/Tutorial 3 - ODEs, Initial Value Problems/J3_ODEIVPsystems.ipynb b/tutorial/julia/Tutorial 3 - ODEs, Initial Value Problems/J3_ODEIVPsystems.ipynb
index d7c73f2..cfa620c 100644
--- a/tutorial/julia/Tutorial 3 - ODEs, Initial Value Problems/J3_ODEIVPsystems.ipynb
+++ b/tutorial/julia/Tutorial 3 - ODEs, Initial Value Problems/J3_ODEIVPsystems.ipynb
@@ -42,9 +42,20 @@
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 1,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "\u001b[32m\u001b[1m Updating\u001b[22m\u001b[39m registry at `C:\\Users\\wilhe\\.julia\\registries\\General`\n",
+ "\u001b[32m\u001b[1m Resolving\u001b[22m\u001b[39m package versions...\n",
+ "\u001b[32m\u001b[1mNo Changes\u001b[22m\u001b[39m to `C:\\Users\\wilhe\\Project.toml`\n",
+ "\u001b[32m\u001b[1mNo Changes\u001b[22m\u001b[39m to `C:\\Users\\wilhe\\Manifest.toml`\n"
+ ]
+ }
+ ],
"source": [
"import Pkg; Pkg.add(\"Plots\"); using Plots"
]
@@ -84,15 +95,123 @@
"\n",
"*Solution:*\n",
"\n",
- "The right hand side function $f(x,y)=x^2y^{1/2}$ for implementing explicit Euler is defined in the function getf1() given below.\n",
+ "The right hand side function $f(x,y)=x^2y^{1/2}$ for implementing explicit Euler is defined in the function f_ex1() given below.\n",
"\n"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 8,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n"
+ ]
+ },
+ "execution_count": 8,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
"# Define the explicit euler algoritm\n",
"function ivp_explicit_euler(f,x0,y0,h,n_out,i_out)\n",
@@ -122,10 +241,8 @@
"x_ex1,y_ex1 = ivp_explicit_euler(f_ex1, x0,y0,h,n_out,i_out)\n",
" \n",
"# Plot y(x) through [0.0,2.0]\n",
- "plot(x_ex1, y_ex1, 'r-o')\n",
- "axis([0.0 2.0 1.0 5.0]) # Set axis limits\n",
- "xlabels!('x') # Set x label to `x`\n",
- "ylabels!('y') # Set y label to `y`"
+ "plot(x_ex1, y_ex1)\n",
+ "xlabel!(\"x\"); ylabel!(\"y\") # Set x label to `x` and y label to `y`"
]
},
{
@@ -154,18 +271,146 @@
"\n",
"\\\\[y_{i+1}^{(k)}:=y_{i+1}^{(k)}-\\frac{R(y_{i+1}^{(k)})}{R'(y_{i+1}^{(k)})}=y_{i+1}^{(k)}-\\frac{y_{i+1}^{(k)} - y_i - h x_{i+1}^2 (y_{i+1}^{(k)})^{1/2}}{1-\\frac{1}{2}h x_{i+1}^2 (y_{i+1}^{(k)})^{-1/2}}\\\\]\n",
"\n",
- "The functions for solving this problem using implicit Euler are defined as getf_IE and getdf_IE in the last section."
+ "The functions for solving this problem using implicit Euler are defined as f_ex1() and df_ex1() in the last section."
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 15,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n"
+ ]
+ },
+ "execution_count": 15,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
"function newton_solve_ie(f, df, y, yold, x, h, tol)\n",
" f_val = f(y,yold,x,h);\n",
- " while abs(f) > tol\n",
+ " while abs(f_val) > tol\n",
" y = y - f_val/df(y,x,h);\n",
" f_val = f(y,yold,x,h);\n",
" end\n",
@@ -181,6 +426,7 @@
" y = y0\n",
" for j = 2:n_out+1\n",
" for i = 1:i_out # Advance i_out steps\n",
+ " yold = y\n",
" y = newton_solve_ie(f, df, y, yold, x, h, tol) # Solve nonlinear system to determine next step\n",
" x = x + h;\n",
" end\n",
@@ -191,17 +437,14 @@
"end\n",
"\n",
"f_ex2(y,yold,x,h) = f = y - yold - h*(x+h)^2*y^(0.5)\n",
- "df_ex2(y,yold,x,h) = y - yold - h*(x+h)^2*y^(0.5\n",
+ "df_ex2(y,x,h) = 1 - 0.5*h*(x+h)^2*y^(-0.5)\n",
"\n",
"tol = 1E-8 # nonlinear solve, absolute convergence tolerance\n",
- "x_ex2,y_ex2 = ivp_implicit_euler(f_ex2,df_ex2,x0,y0,h,n_out,i_out,nl_tol) # run implicit Euler\n",
- "\n",
- "ax = plot(x_ex2, y_ex2, 'b-s') # Make plot for implicit\n",
- "plot!(ax, x_ex1,y_ex1,'r-o') # Make plot for explicit euler \n",
- "axis([0.0 2.0 1.0 6.5]) # Set axis limits\n",
- "xlabels!('x') # Set x label\n",
- "ylabels!('y') # Set y label\n",
- "legend('implicit Euler','explicit Euler','Location','northwest') # Set legend"
+ "x_ex2,y_ex2 = ivp_implicit_euler(f_ex2,df_ex2,x0,y0,h,n_out,i_out,tol) # run implicit Euler\n",
+ "\n",
+ "ax = plot(x_ex2, y_ex2, label = \"Implicit Euler\") # Make plot for implicit solution\n",
+ "plot!(ax, x_ex1, y_ex1, label = \"Explicit Euler\") # Make plot for explicit solution\n",
+ "xlabel!(\"x\"); ylabel!(\"y\") # Set x label to `x` and y label to `y`"
]
},
{
@@ -231,14 +474,150 @@
"\n",
"This particular method is very efficient and easy to program, and represents a good compromise between the number of functional evaluations and the global accuracy. Given its ubiquitousness, this particular fourth-order Runge–Kutta method is normally called **RK4**.\n",
"\n",
- "We consider using RK4 to integrate the function in **Example 1**:"
+ "**Example 3**: We now consider using RK4 to integrate the function in **Example 1**. We'll make use of the right-hand side function previously defined."
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 17,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n"
+ ]
+ },
+ "execution_count": 17,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
"function ivp_RK4(f, x0,y0,h,n_out,i_out)\n",
" xout = zeros(n_out + 1) # Create storage for output\n",
@@ -250,30 +629,25 @@
" for j = 2:n_out+1\n",
" for i = 1:i_out\n",
" k1 = f(x,y);\n",
- " k2 = f(x + 0.5*h,y + 0.5*h*k1);\n",
- " k3 = f(x + 0.5*h,y + 0.5*h*k2);\n",
- " k4 = f(x + h, y + h*k3);\n",
- " y = y + (h/6)*(k1 + 2*k2 + 2*k3 + k4);\n",
- " x = x + h;\n",
+ " k2 = f(x + 0.5*h,y + 0.5*h*k1)\n",
+ " k3 = f(x + 0.5*h,y + 0.5*h*k2)\n",
+ " k4 = f(x + h, y + h*k3)\n",
+ " y = y + (h/6)*(k1 + 2*k2 + 2*k3 + k4)\n",
+ " x = x + h\n",
" end\n",
- " xout(j) = x;\n",
- " yout(j) = y;\n",
+ " xout[j] = x\n",
+ " yout[j] = y\n",
" end\n",
" return xout, yout\n",
"end\n",
"\n",
- "[xout,y_RK4] = ivp_RK4(x0,y0,h,n_out,i_out); % RK4\n",
- "\n",
- "% Plot y(x) through [0.0,2.0]\n",
- "plot(tout,y_RK4,'k-v') % RK4\n",
- "hold on\n",
- "plot(x_IE,y_IE,'b-s') % implicit Euler\n",
- "plot(x_EE,y_EE,'r-o') % explicit Euler\n",
- "axis([0.0 2.0 1.0 6.5]); % Set axis limits\n",
- "xlabel('x') % Set x label\n",
- "ylabel('y') % Set y label\n",
- "legend('RK4','implicit Euler','explicit Euler','Location','northwest') % Set legend\n",
- "hold off "
+ "x_ex3,y_ex3 = ivp_RK4(f_ex1,x0,y0,h,n_out,i_out) # RK4\n",
+ "\n",
+ "# Plot y(x) through [0.0,2.0]\n",
+ "ax = plot(x_ex3, y_ex3, label = \"RK4\") # Make plot for RK4 solution\n",
+ "plot!(ax, x_ex2, y_ex2, label = \"Implicit Euler\") # Make plot for implicit solution\n",
+ "plot!(ax, x_ex1, y_ex1, label = \"Explicit Euler\") # Make plot for explicit solution\n",
+ "xlabel!(\"x\"); ylabel!(\"y\") # Set x label to `x` and y label to `y` "
]
},
{
@@ -281,7 +655,10 @@
"metadata": {},
"source": [
"Again, the truncation errors accumulate through integration resulting in the deviation of the trajactories using different methods. Now, try to reduce the deviation by using smaller step size:\n",
- "Interactive! Reduce the step size to h = 0.01 in the above code, you also need to modify the number of output points n_out for consistency.\n",
+ "\n",
+ "
\n",
+ "INTERACTIVE! Reduce the step size to h = 0.01 in the above code, you also need to modify the number of output points n_out for consistency. How does this affect the behavior of the plots?\n",
+ "
\n",
"\n",
" \n",
"\n",
@@ -343,37 +720,215 @@
"The initial conditions are $y_1(0)=1$ and $y_2(0)=2$. Plot the profiles of $y_1(x)$ and $y_2(x)$ through time horizon $x \\in [0,0.5]$ using both methods.\n",
"\n",
"**Solution:** \n",
- "We first consider to use RK4 for integration: we define the variable in the vector form $\\mathbf y = (y_1,y_2)$, the right hand side function is then given by $\\mathbf f(\\mathbf y)$ and defined as getf_E2 which is written in the last section."
+ "We first consider to use RK4 for integration: we define the variable in the vector form $\\mathbf y = (y_1,y_2)$, the right hand side function is then given by $\\mathbf f(\\mathbf y)$ is defined as f_ex2."
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 19,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "f_ex2 (generic function with 2 methods)"
+ ]
+ },
+ "execution_count": 19,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
- "clear\n",
- "t0 = 0.0; % initial condition\n",
- "y0 = [1.0;2.0]; % initial condition\n",
- "h = 1e-3; % step size\n",
- "n_out = 50; % number of output points\n",
- "i_out = 10; % frequency\n",
- "[tout,y_RK4] = ivp_RK4_sys(t0,y0,h,n_out,i_out); % RK4\n",
- "% getf_E2(y) is required\n",
- "\n",
- "% Plot\n",
- "plot(tout,y_RK4(:,1),'k-o')\n",
- "hold on\n",
- "plot(tout,y_RK4(:,2),'k-x')\n",
- "axis([0.0 0.5 0.0 5.0])\n",
- "xlabel('t')\n",
- "ylabel('y')\n",
- "legend('y_1','y_2','Location','northwest')\n",
- "hold off"
+ "f_ex2(y) = [y[1]*y[2]^2; (y[1]-y[2])^2]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We may then define the RK4 integration scheme (for ODE systems) and use it to integrate the differential equation."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 36,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n"
+ ]
+ },
+ "execution_count": 36,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "function ivp_RK4_sys(f,x0,y0,h,n_out,i_out)\n",
+ " xout = zeros(n_out+1); xout[1] = x0\n",
+ " yout = zeros(n_out+1,length(y0)); yout[1,:] = y0\n",
+ " x = x0; y = y0;\n",
+ " for j = 2:n_out+1\n",
+ " for i = 1:i_out\n",
+ " k1 = f(y)\n",
+ " k2 = f(y + 0.5*h*k1)\n",
+ " k3 = f(y + 0.5*h*k2)\n",
+ " k4 = f(y + h*k3)\n",
+ " y = y + (h/6)*(k1 + 2*k2 + 2*k3 + k4)\n",
+ " x = x + h\n",
+ " end\n",
+ " xout[j] = x\n",
+ " yout[j,:] = y\n",
+ " end\n",
+ " return xout, yout\n",
+ "end\n",
+ " \n",
+ "x0 = 0.0 # initial condition\n",
+ "y0 = [1.0;2.0]\n",
+ "h = 1e-3 # step size\n",
+ "n_out = 50 # number of output points\n",
+ "i_out = 10 # frequency to save point\n",
+ "x_RK4s, y_RK4s = ivp_RK4_sys(f_ex2,x0,y0,h,n_out,i_out) # RK4\n",
+ "\n",
+ "ax = plot(x_RK4s, y_RK4s[:,1], label = \"RK4, y[1]\", # Make plot for component 1 of RK4 solution\n",
+ " xlim = (0.0, 0.5), ylim = (0.0, 5.9)) \n",
+ "plot!(ax, x_RK4s, y_RK4s[:,2], label = \"RK4, y[2]\") # Make plot for component 2 of implicit solution\n",
+ "xlabel!(\"x\"); ylabel!(\"y\") # Set x label to `x` and y label to `y` "
]
},
{
- "attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
@@ -391,36 +946,247 @@
"-2h\\left(y_{1,i+1} -y_{2,i+1} \\right) & 1+2h\\left(y_{1,i+1} -y_{2,i+1} \\right)\n",
"\\end{array}\\right\\rbrack \\\\]\n",
"\n",
- "The residual getR and the Jaconian function getJ are written in the last section"
+ "The residual R_ex2! and the Jaconian function J_ex2! then specified as:"
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 32,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "J_ex2! (generic function with 1 method)"
+ ]
+ },
+ "execution_count": 32,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
- "clear \n",
- "t0 = 0.0; % initial condition\n",
- "y0 = [1.0;2.0]; % initial condition\n",
- "h = 1e-3; % step size\n",
- "n_out = 50; % number of output points\n",
- "i_out = 10; % frequency\n",
- "[x_IE,y_IE] = ivp_implicit_euler_sys(t0,y0,h,n_out,i_out); % implicit Euler\n",
- "% getR(y) and getJ(y) is required\n",
- "\n",
- "plot(x_IE(1:49),y_IE(1:49,1),'bs')\n",
- "hold on\n",
- "plot(x_IE(1:49),y_IE(1:49,2),'b-+')\n",
- "axis([0.0 0.5 0.0 5.0])\n",
- "xlabel('t')\n",
- "ylabel('y')\n",
- "legend('y_1','y_2','Location','northwest')\n",
- "hold off"
+ "function R_ex2!(out, y, yold, h)\n",
+ " out[1] = y[1] - yold[1] - h*y[1]*y[2]^2\n",
+ " out[2] = y[2] - yold[2] - h*(y[1] - y[2])^2\n",
+ " return nothing\n",
+ "end\n",
+ "\n",
+ "function J_ex2!(out, y, h)\n",
+ " out[1,1] = 1 - h*y[2]^2\n",
+ " out[1,2] = -2*h*y[1]*y[2]\n",
+ " out[2,1] = -2*h*(y[1]-y[2])\n",
+ " out[2,2] = 1 + 2*h*(y[1]-y[2])\n",
+ " return nothing\n",
+ "end"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We may then define the implicit Euler integration scheme (for ODE systems) and use it to integrate the differential equation."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 56,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Failed to converge in specified interation limit.\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n"
+ ]
+ },
+ "execution_count": 56,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "using LinearAlgebra: norm\n",
+ "\n",
+ "function ivp_implicit_euler_sys(R!,J!,x0,y0,h,n_out,i_out,tol)\n",
+ " xout = zeros(n_out+1); xout[1] = x0\n",
+ " yout = zeros(n_out+1,length(y0)); yout[1,:] = y0\n",
+ " ny = length(y0)\n",
+ " R = zeros(ny)\n",
+ " J = zeros(ny,ny)\n",
+ " x = x0; y = y0\n",
+ " for j = 2:n_out+1\n",
+ " for i = 1:i_out\n",
+ " yold = y\n",
+ " R!(R,y,yold,h)\n",
+ " count = 1\n",
+ " while norm(R) > tol\n",
+ " J!(J,y,h)\n",
+ " del = -J\\R\n",
+ " y += del\n",
+ " R!(R,y,yold,h)\n",
+ " count += 1\n",
+ " if count > 20\n",
+ " println(\"Failed to converge in specified interation limit.\")\n",
+ " return xout,yout\n",
+ " end\n",
+ " end\n",
+ " x += h\n",
+ " end\n",
+ " xout[j] = x\n",
+ " yout[j,:] = y\n",
+ " end\n",
+ " return xout,yout\n",
+ "end\n",
+ "\n",
+ "x0 = 0.0 # initial condition\n",
+ "y0 = [1.0;2.0]\n",
+ "h = 1e-3 # step size\n",
+ "n_out = 50 # number of output points\n",
+ "i_out = 10 # frequency to save points\n",
+ "tol = 1e-6\n",
+ "x_IEs, y_IEs = ivp_implicit_euler_sys(R_ex2!,J_ex2!,x0,y0,h,n_out,i_out,tol) # implicit Euler\n",
+ "\n",
+ "ax = plot(x_IEs[1:49], y_IEs[1:49,1], label = \"Implicit Euler, y[1]\", # Make plot for component 1 of RK4 solution\n",
+ " xlim = (0.0, 0.5), ylim = (0.0, 5.0)) \n",
+ "plot!(ax, x_IEs[1:49], y_IEs[1:49,2], label = \"Implicit Euler, y[2]\") # Make plot for component 2 of implicit solution\n",
+ "xlabel!(\"x\"); ylabel!(\"y\") # Set x label to `x` and y label to `y`"
]
},
{
- "attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
@@ -460,37 +1226,225 @@
"1\n",
"\\end{array}\\right\\rbrack \\\\]\n",
"\n",
- "and the initial condition is given by $\\mathbf y_0=(1,0,0,0)$. For solving this problem, we consider to use RK4 algorithm."
+ "and the initial condition is given by $\\mathbf y_0=(1,0,0,0)$. For solving this problem, we consider to use RK4 algorithm and define the right-hand side function below."
]
},
{
"cell_type": "code",
- "execution_count": null,
+ "execution_count": 71,
"metadata": {},
- "outputs": [],
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "f_CS_param (generic function with 1 method)"
+ ]
+ },
+ "execution_count": 71,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
"source": [
- "clear\n",
- "F = 100; % Feed flow rate (L/h)\n",
- "v = 100; % Reactor volume (L)\n",
- "CA0 = 1; % Concentration of A in the feed (mol/L)\n",
- "t0 = 0.0; % initial condition\n",
- "y0 = [1.0;0.0;0.0;0.0]; % initial condition\n",
- "h = 1e-3; % step size\n",
- "n_out = 50; % number of output points\n",
- "i_out = 100; % frequency\n",
- "[tout,yout] = ivp_RK4_CS(t0,y0,h,n_out,i_out,F,v,CA0); % RK4\n",
- "% getf_CS(y) is required\n",
- "\n",
- "% Plot\n",
- "plot(tout,yout(:,1),'b-')\n",
- "hold on\n",
- "plot(tout,yout(:,2),'r-')\n",
- "plot(tout,yout(:,3),'k-')\n",
- "axis([0.0 5.0 0.0 1.0])\n",
- "xlabel('t')\n",
- "ylabel('y')\n",
- "legend('C_A','C_B','C_C','Location','northeast')\n",
- "hold off"
+ "function f_CS_param(y,F,v,CA0)\n",
+ " # process specifications and constants\n",
+ " a = 0.1\n",
+ " w = 2*pi # rads/s\n",
+ " k1 = 3 # h^-1\n",
+ " k2 = 1 # h^-1\n",
+ "\n",
+ " # residual\n",
+ " [F/v*(CA0*(a*cos(w*y[4]) + 1) - y[1]) - k1*y[1];\n",
+ " -F/v*y[2] + k1*y[1] - k2*y[2];\n",
+ " -F/v*y[3] + k2*y[2];\n",
+ " 1]\n",
+ "end "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 77,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n"
+ ]
+ },
+ "execution_count": 77,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "F = 100 # Feed flow rate (L/h)\n",
+ "v = 100 # Reactor volume (L)\n",
+ "CA0 = 1 # Concentration of A in the feed (mol/L)\n",
+ "\n",
+ "f_CS(y) = f_CS_param(y,F,v,CA0) # Fix parameters in function\n",
+ "\n",
+ "x0 = 0.0 # initial condition\n",
+ "y0 = [1.0;0.0;0.0;0.0]\n",
+ "h = 1e-3 # step size\n",
+ "n_out = 50 # number of output points\n",
+ "i_out = 100 # save frequency\n",
+ "xout_Cs,yout_Cs = ivp_RK4_sys(f_CS,x0,y0,h,n_out,i_out) # RK4\n",
+ "\n",
+ "ax = plot(xout_Cs, yout_Cs[:,1], label = \"C_A\", # Make plot for component 1 of RK4 solution\n",
+ " xlim = (0.0, 5.0), ylim = (0.0, 1.0)) \n",
+ "plot!(ax, xout_Cs, yout_Cs[:,2], label = \"C_B\") # Make plot for component 2 of implicit solution\n",
+ "plot!(ax, xout_Cs, yout_Cs[:,3], label = \"C_C\") # Make plot for component 2 of implicit solution\n",
+ "xlabel!(\"t\"); ylabel!(\"C\") # Set x label to `x` and y label to `y`"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "
\n",
+ "INTERACTIVE! Now try manipulating the volumetric flowrate of the feed *F*. What flowrate is sufficient to suppress the observed oscillations?\n",
+ "
"
]
},
{
diff --git a/tutorial/matlab/Tutorial 3 - ODEs, Initial Value Problems/M3_ODEIVPsystems.mlx b/tutorial/matlab/Tutorial 3 - ODEs, Initial Value Problems/M3_ODEIVPsystems.mlx
index d191fbf..6680edc 100644
Binary files a/tutorial/matlab/Tutorial 3 - ODEs, Initial Value Problems/M3_ODEIVPsystems.mlx and b/tutorial/matlab/Tutorial 3 - ODEs, Initial Value Problems/M3_ODEIVPsystems.mlx differ
![](\"CSTR_ChloroBenzene.png\")
\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",
+ "\\\\[\n",
+ "\\begin{array}{l}\n",
+ "F_1 -F_2 +F_7 =0 \\\\\n",
+ "F_2 - k_1 y_{3,A}v / (y_{3,A} {\\hat{v} }_A +y_{3,B} {\\hat{v} }_B +y_{3,C} {\\hat{v} }_C) - y_{3,A} F_3 =0\\\\\n",
+ "v(k_2 y_{3,B} -k_1 y_{3,A}) / (y_{3,A} {\\hat{v} }_A +y_{3,B} {\\hat{v} }_B +y_{3,C} {\\hat{v} }_C) + y_{3,B} F_3 =0\\\\\n",
+ "k_2 y_{3,B}v / (y_{3,A} {\\hat{v} }_A +y_{3,B} {\\hat{v} }_B +y_{3,C} {\\hat{v} }_C)-y_{3,C} F_3 =0 \\\\\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",
+ "F_4 -F_5 -F_6 =0\\\\\n",
+ "y_{4,B} F_4 -F_5 =0\\\\\n",
+ "y_{3,A} +y_{3,B} +y_{3,C} -1=0\\\\\n",
+ "y_{4,B} +y_{4,C} -1=0\n",
+ "\\end{array}\n",
+ "\\\\]\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. \n",
+ "\n",
+ "