.. sectionauthor:: Yasuyuki SHIMIZU ############################################## The higher-oder upwind difference scheme ############################################## Interpolation of variable distribution between grids ========================================================= Let us assume that :math:`f(t)` is the distribution of the variable, :math:`f`, of the triangular distribution at time :math:`t`, shown as the solid black line in :numref:`2dcip_1`. and that the state of the advection of the distribution: a velocity of :math:`u(>0)`, and to the direction of :math:`x`. Also, at time :math:`t`, the vertex of triangle, or point :math:`A`, is assumed to be located at the :math:`i` grid point. The distance of triangular distribution moved during the time, :math:`\Delta t`, is :math:`u \Delta t`. Therefore, the distribution supposed to move to the red broken line in the figure and the peak, or vertex, moves to the point :math:`A'` shown by the red circle. .. _2dcip_1: .. figure:: images/01/2dcip_1.png :align: center :width: 70% :: Advection of variable (1) However, actually as shown by :numref:`2dcip_1`, point :math:`A'` is not on the calculation grid, but is in between the :math:`i` th grid and :math:`i+1` th grid. This means that we cannot express the triangular distribution that has moved and whose vertex is point :math:`A'` with discrete values. When we try to express this with a conventional grid point, what we obtain is the shape shown by the blue line in :numref:`2dcip_2`. The peak point :math:`A'` can only be expressed as a shape whose peak has damped, such as point :math:`B_1` or point :math:`B_2`. This is a graphic description of the numerical diffusion that was explained in :ref:`art_dif` of the previous chapter. .. _2dcip_2: .. figure:: images/01/2dcip_2.png :align: center :width: 70% :: Advection of variable (2) The above explanation is for the case where time moves from :math:`t` to :math:`t+\Delta t`. Let us replace this with the case where time moves from :math:`t-\Delta t` to :math:`t`. As shown in :numref:`2dcip_3`, in order for the distribution of :math:`f` at time :math:`t` to be correctly represented as a triangular distribution whose peak is :math:`A`, it is necessary for the triangular distribution at time :math:`t-\Delta t` with a peak of :math:`A'` to be correctly expressed. However, because point :math:`A'` is not on a grid point, it is impossible to express this shape by means of discrete values. This means that the shape of distribution at time :math:`t-\Delta t` must not have discrete values but must have some form of a continuous function. This is the origin of the constrained interpolation profile (CIP) method explained as follows. .. _2dcip_3: .. figure:: images/01/2dcip_3.png :width: 70% :align: center :: Advection of variable (3) .. _label1: The basic concept of the `CIP method `_ =============================================================================================== Let us review the 1D advection equation. .. math:: :label: kiso_11 \cfrac{\partial f}{\partial t}+u\cfrac{\partial f}{\partial x}=0 In :eq:`kiso_1` of the previous chapter, the advection velocity was expressed as :math:`c`, but here we express it as :math:`u`. (There is no special reason for this change of expression. However, because :math:`u, v, w` are often used for expressing flow velocity, here, we use :math:`u`.) The numerical calculation of the 1D advection equation in Equation :eq:`kiso_11` has the same meaning as the estimation of :math:`f(x,t+\Delta t)`. As was explained in the previous section, the normal upwind difference scheme is applied to estimate by using :math:`f(x-\Delta x,t)` and :math:`f(x,t)` As mentioned above, this method provides a stable calculation, but it is not accurate because the calculation results diffuse (numerical diffusion). Here, I introduce the constrained interpolation profile method (`the CIP method `_) developed by Yebe, et al. [Ref:1]_. This higher-order difference scheme replaces the 1D upwind difference scheme. Originally Equation :eq:`kiso_11` means that the distribution profile of variable :math:`f` is conveyed at a speed of :math:`u`. Therefore, the distance of movement of the distribution profile, :math:`f`, during the time, :math:`\Delta t`, should be :math:`u \Delta t`. This means that the future :math:`f=f(x, t+\Delta t)` should be the same as the current :math:`f=f(x-u\Delta t, t)` (Refer to :numref:`trans`). .. _trans: .. figure:: images/01/fig1-2.png :width: 50% :align: center Conveying ::math:`f(x,t+\Delta t)` Instead of estimating the future :math:`f(x,t+\Delta t)`, if we can estimate the current :math:`f(x-u\Delta t, t)`, that value can be used as an estimate of the future value. This can be expressed by the following formula. .. math:: :label: yosoku f(x,t+\Delta t) \simeq f(x-u\Delta t, t) Method for estimating the variable distribution form between grids ===================================================================== Introduction of a 3rd order function ----------------------------------------------------------- There are many methods for estimating :math:`f(x-u\Delta t, t)`. The issue is that because :math:`f` is given as discrete values, the known quantity of :math:`f` is the point, ● of :numref:`trans`. From the two ● points, :math:`f(x-u\Delta t, t)` of ○ needs to be estimated. The simplest method is linking the two dots by a line and proportionally allocating the quantity. However, when :math:`f` rapidly changes or when it increases and decreases between the two points, it is natural to express the estimated value by using a 3rd order function. .. _moshiki: .. figure:: images/01/fig1-3.png :width: 50% :align: center : Schematic diagram of the CIP method As shown by :numref:`moshiki` discrete values are expressed in the direction of :math:`x` as :math:`i` and in the direction of :math:`t` as :math:`n`. We denote :math:`f` at points A and B as :math:`f(i-1,n)` and :math:`f(i,n)`. Assuming that the distribution of :math:`f` between A and B is expressed as a 3rd order function and that the distance from B to A is expressed as :math:`X`, then we have... .. math:: :label: fvalue F=a_1 X^3 + a_2 X^2 + a_3 X + a_4 This is a 3rd order function that has four parameters, :math:`a_1 \sim a_4`. Here, note that :math:`f` gives discrete values defined for each grid point, while :math:`F` is a 3D continuous function that interpolates :math:`f` of each grid point between the grids. To determine the four parameters in Equation :eq:`fvalue`, the following four conditions are necessary. Two of them are given by... .. math:: :label: bc0 \mbox{at point B:} & F_{(X=0)}=f(i,n) \\ \mbox{at point A:} & F_{(X=\Delta x)} = f(i-1,n) What about the remaining two conditions? We can think about them in many ways (for example, by using :math:`f(i-2,n)` or :math:`f(i+1,n)`). In the CIP method, we use the gradient of :math:`f` in the direction of the :math:`x` axis. Thus, .. math:: :label: bc1 \mbox{at B:} & \displaystyle{{\partial F}\over{\partial x}}_{(X=0)}= {{\partial f}\over{\partial x}} (i,n) \\ \mbox{at A:} & \displaystyle{{\partial F}\over{\partial x}}_{(X=\Delta x)} = {{\partial f}\over{\partial x}}(i-1,n) By using equations :eq:`bc0` :math:`\sim`:eq:`bc1`, only with information given by :math:`i` and :math:`i-1`, :math:`F`, the distribution of :math:`f` between A and B can be estimated. This is one of the major advantages of the CIP method. Equation :eq:`fvalue` and its partial differentiation in terms of :math:`x` are... .. math:: :label: fprof F_x= {{\partial F}\over{\partial x}} =3 a_1 X^2 + 2 a_2 X +a_3 For the above equation, by applying the four conditions of equations :eq:`bc0` :math:`\sim` :eq:`bc1`, we can determine the parameters, :math:`a_1 \sim a_4`, and obtain :math:`F`. However, it seems a little complicated. Then, among equations :eq:`bc0` :math:`\sim` :eq:`bc1`, the conditions of point B (:math:`X=0`) automatically satisfied by .. math:: :label: fxnew F(X)=\left[ (aX + b)X + f_x(i,n) \right] X + f(i,n) Now, let us use this equation. Note that .. math:: :label: fnew1 f_x(i,n) = \displaystyle{{\partial f}\over{\partial x}} (i,n) Convert Equation :eq:`fxnew` to .. math:: :label: fnew2 F(X)=aX^3+bX^2+f_x(i,n)X+f(i,n) And when we perform the partial differential in terms of :math:`x`, we have... .. math:: :label: fnew3 F_x(X)=3aX^2+2bX+f_x(i,n) Equation :eq:`fnew2` and Equation :eq:`fnew3` automatically satisfy the conditions of point B (:math:`X=0`) of Equation :eq:`bc0` and Equation :eq:`bc1`. Now, let us apply the conditions of point B (:math:`X=\Delta x`) of Equation :eq:`bc0` and Equation :eq:`bc1` to Equation :eq:`fnew2` and Equation :eq:`fnew3`. Then, we obtain... .. math:: :label: ab1 a\Delta x^3 + b \Delta x^2 + f_x(i,n) \Delta x + f(i,n) = f(i-1,n) \\ 3a\Delta x^2 + 2b\Delta x + f_x(i,n) = f_x(i-1,n) Therefore, we obtain... .. math:: :label: ab2 a=\displaystyle{{ \left[(f_x(i-1,n)+f_x(i,n) \right] \Delta x + 2 \left[f(i,n)-f(i-1,n) \right] }\over{\Delta x^3}} \\ b=\displaystyle{{ 3\left[(f(i-1,n)-f(i,n) \right]- \left[f_x(i-1,n)+2f_x(i,n) \right] \Delta x }\over{\Delta x^2}} Problem (1) ^^^^^^^^^^^^^^ .. list-table:: * - Drive Equation :eq:`ab2`. Summary of a method for estimating the distribution equation of link variables by using the CIP method. -------------------------------------------------------------------------------------------------------------- In summary, the procedure for applying the CIP method to Equation :eq:`kiso_11` in order to obtain :math:`f(x,t+\Delta t)` is as follows. 1) Obtain :math:`a, b` for Equation :eq:`ab2`. 2) Substitute :math:`a, b` into Equation :eq:`fxnew` to transform it to :math:`X=u\Delta t` so as to obtain :math:`F(X=u\Delta t)=f(x-u\Delta t, t)`. 3) At the same time, use Equation :eq:`fnew3` to obtain :math:`F_x(X=u\Delta t)=f_x(x-u\Delta t, t)`. 4) Obtain :math:`f(x,t+\Delta t)` for Equation :eq:`yosoku`. Execute this procedure for :math:`i=1 \sim N_i`, and when it finishes, update the time. The Python code and calculation examples ---------------------------------------------- An example Python code for the numerical solution of the 1D advection equation by the CIP method is given here. https://github.com/YasuShimizu/Python-Code/blob/main/cip.py Performing this calculation yields the following results. Where, :math:`x_b` = 10, :math:`x_p` = 10, :math:`f_p` = 0.5, :math:`u` = 0.5． .. _cip_tri: .. figure:: images/01/cip.gif :width: 90% :align: center : The 1D advection equation (CIP method) The CIP method can thus be used to calculate how the triangular distribution in the initial conditions is not disrupted, and how it moves downstream at a constant speed. The results are almost identical to the theoretical solution (:numref:`adv_anim`) in the Appendix. Advection equations with varying advection directions ============================================================== Whereas :ref:`label1` above described a calculation method for the case of constant advection direction, this section describes a method for the case where the advection direction varies. In the 1D advection equation, :eq:`kiso_11`, the advection velocity, :math:`u`, was treated as constant in the previous explanations. Here, :math:`u` is not a constant, but is a quantity that varies temporally and spatially. In general, in fluid, :math:`u` is positive or negative depending on the direction of flow. The algorithm must be able to cope with both :math:`u>0` and :math:`u<0`. When *u<0*, ------------------ .. _moshiki_u-: .. figure:: images/01/fig1-4.png :width: 50% :align: center : Schematic for when :math:`u<0`. As shown in the above figure, when :math:`u<0`, we assume that... .. math:: :label: X_nega X=-u\Delta t and that the distribution of :math:`f` between :math:`i` and :math:`i+1` is expressed with the following equation. Where, :math:`X>0`. .. math:: :label: X_nega_1 F(X)=\left[ (aX + b)X + f_x(i) \right] X + f(i) Here, the differentiation of :math:`F(X)` is .. math:: :label: X_nega_2 F'(X)=3aX^2+2bX+f_x(i) Where, :math:`f_x` is .. math:: :label: X_nega_3 f_x = {{\partial f}\over{\partial x}} :math:`F` and :math:`F'` derive from the conditions of... .. math:: :label: X_nega_4 X=0 \mbox{: } & F(X)=f(i) \mbox{ and } F'(X)=f_x(i) \\ X=\Delta x \mbox{: } & F(X)=f(i+1) \mbox{ and } F'(X)=f'(i+1) as follows. .. math:: :label: X_nega_5 a=\displaystyle{{ \left[(f_x(i+1)+f_x(i) \right] \Delta x + 2 \left[f(i)-f(i+1) \right] }\over{\Delta x^3}} \\ b=\displaystyle{{ 3\left[(f(i+1)-f(i) \right]- \left[f_x(i+1)+2f_x(i) \right] \Delta x }\over{\Delta x^2}} Problem (2) ^^^^^^^^^^^^^^^^ .. list-table:: * - Drive Formula :eq:`X_nega_5`. When *u>0*, ----------------------------------------------- .. figure:: images/01/fig1-5.png :width: 50% :align: center : Schematic diagram for when :math:`u>0` As shown in the above figure, when :math:`u>0`, we assume that… .. math:: :label: X_posi_1 X = - u \Delta t Then, despite the assumption that :math:`X = u \Delta t (>0)` in the previous section, here, we assume u such that :math:`X<0`, as in :eq:`X_posi_1`. The reasons will be explained later, but the assumption enables us to use the same algorithm regardless of whether :math:`u<0` or :math:`u>0`. For :math:`u>0`, the distribution of :math:`f` between :math:`i` and :math:`i-1` is expressed by the following equation. This is completely the same form as Equation :eq:`X_nega_1`. .. math:: :label: X_posi_2 F(X)=\left[ (aX + b)X + f_x(i) \right] X + f(i) Then, the differentiation of :math:`F(X)` is expressed as follows. .. math:: :label: X_posi_3 F'(X)=3aX^2+2bX+f_x(i) From the following conditions, .. math:: :label: X_posi_4 X=0 \mbox{: } & F(X)=f(i) \mbox{ and } F'(X)=f_x(i) \\ X= -\Delta x \mbox{t: } & F(X)=f(i-1) \mbox{ and } F'(X)=f'(i-1). :math:`F` and :math:`F'` give the following formulas. .. math:: :label: X_posi_5 a=\displaystyle{{ \left[(f_x(i-1)+f_x(i) \right] \Delta x - 2 \left[f(i)-f(i-1) \right] }\over{\Delta x^3}} \\ b=\displaystyle{{ 3\left[(f(i-1)-f(i) \right]+ \left[f_x(i-1)+2f_x(i) \right] \Delta x }\over{\Delta x^2}} Problem (3) ^^^^^^^^^^^^^^^ .. list-table:: * - Drive Formula :eq:`X_posi_5`. How can we make it usable for both *u>0* and *u<0*? -------------------------------------------------------- Formula :eq:`X_nega_5` is when :math:`u<0`, and Formula :eq:`X_posi_5` is when :math:`u>0`. These pairs are very similar. Then, we use the sign function (also known as the signal function or the signum function) that is defined as follows. .. math:: :label: sig \mbox {sign}(Y) & =1 & ,Y>0 \\ \mbox {sign}(Y) & =-1 & ,Y<0 In addition, we define :math:`i_m` as follows. .. math:: :label: sig_1 i_m = i - \mbox{sign}(u) Then, we have... .. math:: :label: sig_2 u>0 \mbox{: } i_m & = i-1 \\ u<0 \mbox{: } i_m & = i+1 Therefore, by assuming the following equations, instead of Equation :eq:`X_nega_5` and Equation :eq:`X_posi_5`... .. math:: :label: ab a=\displaystyle{{ \left[(f_x(i_m)+f_x(i) \right] \Delta x \cdot \mbox{sign}(u)- 2 \left[f(i)-f(i_m) \right] }\over{\Delta x^3 \cdot \mbox{sign}(u)}}\\ b=\displaystyle{{ 3\left[(f(i_m)-f(i) \right]+ \left[f_x(i_m)+2f_x(i) \right] \Delta x \cdot \mbox{sign}(u)}\over{\Delta x^2}} Then, we can use the same equations as for when :math:`u<0` and for when :math:`u>0`. Example ^^^^^^^^^ .. list-table:: * - Use the CIP method to calculate the 1D advection equation for the case where the advection velocity :math:`u` exhibits simple harmonic motion in the positive or negative direction, expressed in the following equation. * - The initial conditions are the same as those for :numref:`triangle`. .. math:: :label: uprof u=u_0 \sin \left( \cfrac{2\pi}{T_l} t \right) Where, :math:`u_0` is the maximum value of the oscillatory flow velocity, and :math:`T_l` is the period of the oscillatory flow velocity. The Python code and calculation examples ---------------------------------------------- When the advection direction of the advection equation oscillates positively and negatively, an example of the Python code for numerical calculation by means of the CIP method is below. https://github.com/YasuShimizu/Python-Code/blob/main/cip_pm.py With the code, assuming that :math:`u_0` = 2.0, :math:`T_l` = 100, :math:`x_b` = 10, :math:`x_p` = 20, and :math:`f_p` = 0.5, then we have :numref:`pm`, which shows the calculation result. .. _pm: .. figure:: images/01/cip_pm.gif :width: 50% :align: center : Calculation results when the advection direction oscillates positively and negatively An advection equation with the source term ================================================ This section describes the splitting method, which is indispensable when performing numerical calculations for fluids. Fluid equations rarely have zero on the right side, as in the Equation :eq:`kiso_11`, and many have a value on the right side. This is generally expressed as follows. .. math:: :label: sr1 \cfrac{\partial f}{\partial t}+u\cfrac{\partial f}{\partial x}=G For the equations of motion of a fluid, :math:`G` includes various external forces (or source terms), such as the surface slope, the bed slope, the friction, etc. in short, it is other elements than the advection term. For convenience, the Equation :eq:`sr1` is divided into the following two equations. .. math:: :label: bunri1 {{\partial f}\over{\partial t}} = G .. math:: :label: bunri2 {{\partial f}\over{\partial t}}+u{{\partial f}\over{\partial x}} = 0 :math:`f(t)` and :math:`f(t+\Delta t)` are set as :math:`f^n` and :math:`f^{n+1}`, respectively. The calculation of :math:`f^n \rightarrow f^{n+1}` by Equation :eq:`sr1` is performed as follows. First, we use Equation :eq:`bunri1` to calculate :math:`f^n \rightarrow \widetilde{f}`. Second, calculate :math:`\widetilde{f} \rightarrow f^{n+1}` with Equation :eq:`bunri2`. This two-phase calculation is called the splitting method. Only after splitting the calculation like this, can the CIP method be applied in the form of Equation :eq:`bunri2`. These two phases are called the non-advection phase and the advection phase. In the CIP method, in addition to the time variation of :math:`f` the spatial differentiation amount, :math:`f_x`, represented by Equation :eq:`ab` is always required. Therefore, we calculate the differentiation amounts at the same time, as follows. For the non-advection phase, :math:`\displaystyle{{{\partial f}\over{\partial x}}^n \rightarrow \widetilde{{\partial f}\over{\partial x}}}`, For the advection phase, :math:`\displaystyle{\widetilde{{\partial f}\over{\partial x}} \rightarrow {{\partial f}\over{\partial x}}^{n+1}}` In this way, they are calculated separately in two phases. The above is summarized below. .. _bunri_hyo: .. list-table:: Splitting method :header-rows: 1 * - Phase name - Equation used - Updating :math:`f` - Updating :math:`\cfrac{\partial f}{\partial x}` * - **[A]** \ non-advection phase - :math:`\displaystyle{{\partial f}\over{\partial t}}=G` - **(a)** \ :math:`f^n \rightarrow \widetilde{f}` - **(b)** \ :math:`\displaystyle{{{\partial f}\over{\partial x}}^n \rightarrow \widetilde{{\partial f}\over{\partial x}}}` * - **[B]** \ advection phase - :math:`\displaystyle{{{\partial f}\over{\partial t}}+u{{\partial f}\over{\partial x}}}=0` - **(c)** \ :math:`\widetilde{f} \rightarrow f^{n+1}` - **(d)** \ :math:`\displaystyle{\widetilde{{\partial f}\over{\partial x}} \rightarrow {{\partial f}\over{\partial x}}^{n+1}}` A method for using the split-equations to calculation each phase is given below. In the following explanation, if it is necessary to show differences other than the advection term, then the central difference scheme with the highest possible accuracy should be used. Non-advection phase ------------------------- From Equation :eq:`bunri1`, .. math:: :label: na1 {{\widetilde{f}-f^n}\over{\Delta t}} = G^n .. math:: :label: na2 \widetilde{f} = f^n +G^n \Delta t Assuming that the calculation point is :math:`i`, then at... .. math:: :label: na3 \widetilde{f}(i) = f^n(i) +G^n(i) \Delta t :math:`f` of the non-advection phase is updated (:math:`f^n \rightarrow \widetilde{f}`). We differentiate both sides of :eq:`na2` with respect to :math:`x`, then we have... .. math:: :label: na4 \widetilde{{\partial f}\over{\partial x}} ={{\partial f}\over{\partial x}}^n +{{\partial G}\over{\partial x}}^n \Delta t Expressing :math:`\displaystyle{{{\partial G}\over{\partial x}}^n}` in the above equation with the central difference scheme give us... .. math:: :label: na5 \widetilde{{\partial f}\over{\partial x}} ={{\partial f}\over{\partial x}}^n +{{G^n(i+1)-G^n(i-1)}\over{2\Delta x}} \Delta t Then, because of Equation :eq:`na3`, we have... .. math:: :label: na6 G^n(i+1)={{\widetilde{f}(i+1)-f^n(i+1)}\over{\Delta t}} \\ G(i-1)={{\widetilde{f}(i-1)-f^n(i-1)}\over{\Delta t}} Here, let us substitute :eq:`na6` into :eq:`na5`. Then, we have... .. math:: :label: na8 \widetilde{{\partial f}\over{\partial x}}(i)={{\partial f}\over{\partial x}}^n(i)+{{1}\over{2\Delta x}} \left[\widetilde{f}(i+1)-f^n(i+1)-\widetilde{f}(i-1)+f^n(i-1)\right] Therefore, updates of variables, :math:`f`, and their spatial differentiation amount, :math:`\displaystyle{{{\partial f}\over{\partial x}}}`, in the non-advection phase can be obtained as follows. (a) in Table :numref:`bunri_hyo` is given by Equation :eq:`na3`, and (b) in Table :numref:`bunri_hyo` is given by Equation :eq:`na8`. Advection phase --------------------- Updating :math:`f` is performed in (c) of Table :numref:`bunri_hyo`. The CIP methods described so far can be applied without modification. That means that we apply Equation :eq:`X_nega_1` or Equation :eq:`X_posi_2` (they are the same), then we have... .. math:: :label: ir1 {f(i)}^{n+1}=F(X)= \left[(aX+b)X+ \widetilde{{\partial f}\over{\partial x}}(i) \right] X +\widetilde{f}(i) Where, :math:`a` and :math:`b` are given by Equation :eq:`ab`. Next is a method for updating the spatial differentiation amount, :math:`\displaystyle{{\partial f}\over{\partial x}}`, of the advection phase [(d) in the table :numref:`bunri_hyo`]. Perform partial differentiation on both sides of Equation :eq:`bunri2` with respect to :math:`x`. .. math:: :label: ad11 {{\partial \over \partial x}}\left( {{\partial f}\over{\partial t}} \right) + {{\partial \over \partial x}} \left( u {{\partial f}\over{\partial x}} \right) =0 Where, we assume that :math:`\displaystyle{{\partial f}\over{\partial x}}` is :math:`f'`. And develop the second term. Then, we have .. math:: :label: ad12 {{\partial f'}\over{\partial t}} + u{{\partial f'}\over{\partial x}} = -{{\partial u}\over{\partial x}} f' What we are doing here is updating :math:`\displaystyle{{\partial f}\over{\partial x}}` in the advection phase. What we need to do is to use this equation to simply perform :math:`\widetilde{f'} \rightarrow {f'}^{n+1} \left( \displaystyle{ \widetilde{{\partial f}\over{\partial x}} \rightarrow {{\partial f}\over{\partial x}}^{n+1}} \right)`. Now, let us compare Equation :eq:`ad11` and Equation :eq:`sr1`. Equation :eq:`ad12` is equivalent to Equation :eq:`sr1` when substituting its :math:`f` with :math:`f'` and :math:`G` with :math:`\displaystyle{-{{\partial u}\over{\partial x}} f'}`. As we did for Equation :eq:`sr1`, we apply the splitting method to Equation :eq:`ad11`. Therefore, we split the :eq:`ad11` into... .. math:: :label: br_a1 {{\partial f'}\over{\partial t}} + u{{\partial f'}\over{\partial x}} = 0 .. math:: :label: br_a2 {{\partial f'}\over{\partial t}} = -{{\partial u}\over{\partial x}} f' Assuming a middle value, :math:`\widehat{f'}`,　between :math:`\widetilde{f'}` and :math:`{f'}^{n+1}`, we calculate :math:`\widetilde{f'} \rightarrow \widehat{f'}` with Equation :eq:`br_a1`, and :math:`\widehat{f'} \rightarrow {f'}^{n+1}` with Equation :eq:`br_a2`. One is to calculate Equation :eq:`br_a1` by using the CIP method again. However, it seems better to use either Equation :eq:`X_nega_2` (the spatial differentiation amount of :math:`F(X)` ) or Equation :eq:`X_posi_3` (also the spatial differentiation amount of :math:`F(X)` ), which are used in the CIP method. (Equation :eq:`X_nega_2` and Equation :eq:`X_posi_3` have the same form.) Therefore, assuming that... .. math:: :label: br_a3 \widehat{f'}(i) = \widetilde{F'}(X) and calculate with .. math:: :label: br_a4 \widehat{f'}(i) = 3aX^2 + 2bX + \widetilde{f'}(i) Then, we apply the central difference scheme for Equation :eq:`br_a2`, and we have... .. math:: :label: br_a5 {{{f'}^{n+1}(i) - \widehat{f'}(i)}\over{\Delta t}} = - {{u(i+1)-u(i-1))}\over{2\Delta x}} \widehat{f'}(i) Consequently, we calculate using the following formula. .. math:: :label: br_a6 {f'}^{n+1}(i) = \widehat{f'}(i) -{{u(i+1)-u(i-1))}\over{2\Delta x}} \widehat{f'}(i) \Delta t Problem (4) ^^^^^^^^^^^^^^ .. list-table:: * - When :math:`G` on the right side of Equation :eq:`sr1` is :math:`\displaystyle{D{{\partial^2 f}\over{\partial x^2}}}`, .. math:: \displaystyle{{{\partial f}\over{\partial t}}+u{{\partial f}\over{\partial x}}= D{{\partial^2 f}\over{\partial x^2}}} .. list-table:: * - it becomes an advection and diffusion equation. Describe the method to calculate this with the splitting method. The Python code and calculation examples ---------------------------------------------- You can find an example of the Python code for calculating advection-diffusion equations at... https://github.com/YasuShimizu/Python-Code/blob/main/cip_pm_diffusion.py With the code, now assuming that the diffusion coefficient :math:`D=0.5` = 0.5, :math:`u_0` = 2.0, and :math:`T_l` = 100, then the calculation results is shown in :numref:`pm_dif`. It shows how the initially triangular distribution gradually becomes more diffuse as it moves from left to right, and vice versa. .. _pm_dif: .. figure:: images/01/cip_pm_dif.gif :width: 50% :align: center : Advection-diffusion equation calculation results when the advection velocity oscillates positively and negatively A number of high-precision difference methods have been proposed for advection terms, in addition to the CIP method, some of which are presented in the Appendix (:ref:`others`).