Skip to content

Commit

Permalink
Finalize julia tutorial 4
Browse files Browse the repository at this point in the history
  • Loading branch information
@mewilhel
mewilhel committed Dec 21, 2020
1 parent eb19a24 commit 5566a27
Show file tree
Hide file tree
Showing 6 changed files with 297 additions and 511 deletions.
177 changes: 162 additions & 15 deletions tutorial/julia/Tutorial 3 - ODEs, Initial Value Problems/J3_ODEIVPsystems.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@
"metadata": {},
"outputs": [],
"source": [
"import Pkg; Pkd.add(\"Plots.jl\"); using Plots"
"import Pkg; Pkg.add(\"Plots\"); using Plots"
]
},
{
Expand Down Expand Up @@ -293,38 +293,64 @@
"\n",
"with the initial condition $\\mathbf y (0)=\\mathbf y_0$. Here, we see that $\\mathbf f:\\mathbb R^n\\to\\mathbb R^n$ and $\\mathbf y\\in\\mathbb R^n$, thus we have a system of equations that must be solved simultaneously. \n",
" \n",
" Explicit Methods\n",
"## Explicit Methods\n",
" \n",
"The methods that we have described for single equations can be readily adapted for systems of equations, and the accuracy properties of the methods are unchanged. In particular, the adaptation for explicit methods is really straightforward. All that changes is that you have to step through a system of equations instead of a single equation.\n",
"Explicit Euler\n",
" \n",
"### Explicit Euler\n",
"The explicit Euler method in this case has the form\n",
"\n",
"Predictor-corrector\n",
"\\\\[\\mathbf y_{i+1} := \\mathbf y_i + h \\mathbf h (\\mathbf y_i)\\\\]\n",
" \n",
"### Predictor-corrector\n",
"The predictor–corrector method is given by\n",
"\n",
"RK4\n",
"\\\\[\\mathbf y_{i+1}^P:=\\mathbf y_i+h \\mathbf f(\\mathbf y_i)\\\\\n",
"\\mathbf y_{i+1}^C:=\\mathbf y_i+h \\mathbf f(\\mathbf y_{i+1}^P)\\\\\n",
"\\mathbf y_{i+1}:=\\frac{\\mathbf y_{i+1}^P+\\mathbf y_{i+1}^C}{2} \\\\]\n",
" \n",
"### RK4\n",
"The RK4 scheme becomes\n",
"\n",
"\\\\[ \\mathbf k_1:=\\mathbf f_1(\\mathbf y_i)\\\\\n",
"\\mathbf k_2:=\\mathbf f_2(\\mathbf y_i+\\frac{h}{2}\\mathbf k_1)\\\\\n",
"\\mathbf k_3:=\\mathbf f_3(\\mathbf y_i+\\frac{h}{2}\\mathbf k_2)\\\\\n",
"\\mathbf k_4:=\\mathbf f_4(\\mathbf y_i+h \\mathbf k_3)\\\\\n",
"\\mathbf y_{i+1}:=y_i+\\frac{h}{6}(\\mathbf k_1+2\\mathbf k_2+2\\mathbf k_3+\\mathbf k_4) \\\\]\n",
" \n",
"Implicit Euler\n",
"### Implicit Euler\n",
"The concepts for solving systems of ODEs are similar for implicit methods, but a bit more complicated. The form of the implicit Euler method becomes \n",
",\n",
"\n",
"\\\\[\\mathbf y_{i+1} := \\mathbf y_i + h \\mathbf f (\\mathbf y_{i+1})\\\\],\n",
" \n",
"which is a system of nonlinear algebraic equations requiring Newton–Raphson for solving (at each step). We therefore form the residual vector\n",
".\n",
"In the algorithm, we have an outer loop that is updating the values of . At each step , we have an inner loop (with counter ) that involves first solving\n",
"\n",
"and then updating the value for ,\n",
"\\\\[\\mathbf R(y_{i+1})=\\mathbf y_{i+1} - \\mathbf y_i - h \\mathbf f (\\mathbf y _{i+1})=\\mathbf 0\\\\].\n",
"\n",
"In the algorithm, we have an outer loop that is updating the values of $\\mathbf y_i$. At each step $i$, we have an inner loop (with counter $k$) that involves first solving\n",
"\n",
"Example 2: Use RK4 and implicit Euler to solve the system of equations\n",
"\\\\[\\mathbf J ( \\mathbf y_{i+1}^{(k)}) \\mathbf \\delta ^ {(k+1)}=- \\mathbf R ( \\mathbf y_{i+1}^{(k)}) \\\\]\n",
" \n",
"and then updating the value for $\\mathbf y_{i+1}$,\n",
"\n",
"\\\\[\\mathbf y_{i+1}^{(k+1)} := \\mathbf y_{i+1}^{(k)} + \\mathbf \\delta ^ {(k+1)}\\\\]\n",
" \n",
"**Example 2:** Use RK4 and implicit Euler to solve the system of equations\n",
"\n",
"\\\\[\\frac{dy_1}{dt}=y_1y_2^2 \\\\\n",
"\\frac{dy_2}{dt}=(y_1-y_2)^2\\\\]\n",
" \n",
"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",
"The initial conditions are and . Plot the profiles of and through time horizon using both methods.\n",
"Solution: \n",
"We first consider to use RK4 for integration: we define the variable in the vector form , the right hand side function is then given by and defined as getf_E2 which is written in the last section."
"**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."
]
},
{
"cell_type": "markdown",
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"clear\n",
"t0 = 0.0; % initial condition\n",
Expand All @@ -346,6 +372,127 @@
"hold off"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Then, we consider to use implict Euler for integration, the residual can be easily obtained as\n",
"\n",
"\\\\[ \\mathbf{R}=\\left\\lbrack \\begin{array}{c}\n",
"y_{1,i+1} -y_{1,i} -{hy}_{1,i+1} y_{2,i+1}^2 \\\\\n",
"y_{2,i+1} -y_{2,i} -{h\\left(y_{1,i+1} -y_{2,i+1} \\right)}^2 \n",
"\\end{array}\\right\\rbrack = \\mathbf{0} \\\\]\n",
"\n",
"The Jacobian is\n",
"\n",
"\\\\[ \\mathbf{J}=\\left\\lbrack \\begin{array}{cc}\n",
"1-{hy}_{2,i+1}^2 & -2{hy}_{1,i+1} y_{2,i+1} \\\\\n",
"-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"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"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"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice that the solution is changing very rapidly, and this system diverges around . As a result, it is good that we put in a check for the convergence of Newton–Raphson into our program; otherwise, the solution gets stuck in an infinite loop near the point of divergence.\n",
" \n",
"<hr style=\"border:6px solid black\"> </hr>\n",
"\n",
"# <span style=\"color:darkblue\"> Case Study: Series Reactions in a Continuously Stirred Tank Reactor </span>\n",
"\n",
"<img src=\"CSTR_pic.png\" width=\"300\">\n",
"\n",
"As indicated by the above flowchart, species A has a series reaction in a continuous stirred tank reactor (CSTR):\n",
"\n",
"\\\\[\\\\]\n",
"\n",
"$F$ is the volumetric flow rate of the feed stream with pure $A$ ($F=100$ L/h); $v$ is the volume of the reactor ($v=100$ L). The reaction rates for each species are given by:\n",
"\n",
"\\\\[\\frac{dC_A}{dt}=\\frac{F}{v}(C_{A,in}-C_A)-k_1C_A \\\\\n",
"\\frac{dC_B}{dt}=-\\frac{F}{v} C_B + k_1 C_A - k_2 C_B \\\\\n",
"\\frac{dC_C}{dt} = -\\frac{F}{v}C_C + k_2C_B\\\\]\n",
"\n",
"Real processes often have fluctuating feed conditions, which means the feed stream may have a sinusoidal relationship $C_{A,in}=C_{A,0}(A\\cos(\\omega t) + 1)$, the system equations becomes\n",
"\n",
"\\\\[\\frac{dC_A}{dt}=\\frac{F}{v}(C_{A,0}(\\alpha \\cos(\\omega t)+1)-C_A)-k_1C_A \\\\\n",
"\\frac{dC_B}{dt}=-\\frac{F}{v} C_B + k_1 C_A - k_2 C_B \\\\\n",
"\\frac{dC_C}{dt} = -\\frac{F}{v}C_C + k_2C_B\\\\]\n",
"\n",
"The rate constants are $k_1=3\\rm h^{-1}$ and $k_2=1 \\rm h^{-1}$, the amplitude for fluctuations in feed is $A=0.1$, and the angular frequency is given by $\\omega=2\\pi \\frac{\\rm rad}{\\rm s}$. In addition, the initial conditions are given by $C_A(t=0)=C_{A,0}=1 \\frac{\\rm mol}{\\rm L}$, $C_B(t=0)=0 \\frac{\\rm mol}{\\rm L}$ and $C_C(t=0)=0 \\frac{\\rm mol}{\\rm L}$. We want to plot the profile of these species through time $[0.0,5.0]$.\n",
"\n",
"### Solution:\n",
"Immediately, we can see that the system is non-autonomous because of the dependence on $t$ in the first equation. Our first step is to convert it into an autonomous system. Our variable vector is initially $\\mathbf y = (C_A,C_B,C_C)$, by introducing $y_4=t$, the variable vector becomes $\\mathbf y = (C_A,C_B,C_C,t)$. Thus the right hand side function can be rewritten as:\n",
"\n",
"\\\\[\\mathbf{f}\\left(\\mathbf{y}\\right)=\\left\\lbrack \\begin{array}{c}\n",
"\\frac{F}{v}\\left(C_{A,0} \\left(\\alpha \\mathrm{cos}\\left(\\omega y_4 \\right)+1\\right)-y_1 \\right)-k_1 y_1 \\\\\n",
"-\\frac{F}{v}y_2 +k_1 y_1 -k_2 y_2 \\\\\n",
"-\\frac{F}{v}y_3 +k_2 y_2 \\\\\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"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"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down
Loading

0 comments on commit 5566a27

Please sign in to comment.