Linear 1D Advection Equation **************************** Introduction ==================== 1D linear advection equation (so called wave equation) is one of the simplest equations in mathematics. The equation is described as: .. math:: \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 :label: adv_1d where :math:`u(x,t), x \in \mathbb{R}` is a scalar (wave), advected by a nonezero constant :math:`c` during time :math:`t`. The sign of :math:`c` characterise the direction of wave propagation. If :math:`c>0`, the wave propagates in the positive direction of :math:`x`-axis. On the other hand if :math:`c<0` then the wave propagates to the negative direction of `x`-axis. The mangitude of :math:`c` tells us how fast the wave propagates. The exact solution of this equation is given by its initial value. Assuming that the initial value for equation :eq:`adv_1d` is given as :math:`u_0(x) = u(x,0)`, the exact solution of equation :eq:`adv_1d` is: .. math:: u(x,t) = u_{0}(x-ct) :label: ex_adv_1d Numerical Techniques ==================== Although this equation is simple to solve, it can be very useful for learning numerical techniques. We start by introducing the numerical methods. Forward-Time Central-Space (FTCS) Method ---------------------------------------- FTCS is based on central spatial difference scheme and the temporal forward Euler method. Assume that :math:`t` and :math:`x` are descritized uniformly as: .. math:: :label: ti t_{n}= t_{0}+ \Delta t \qquad n=0,...,T x_{i}= x_{0}+ \Delta x \qquad i=0,...,N Let us assume that :math:`u_{i}^{n}:=u(x_i,t_n)`. Then using applying central spatial difference and the temporal forward Euler methods to equation :eq:`adv_1d`, we have: .. math:: \frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t} + c \frac{u_{i-1}^{n}-u_{i+1}^{n}}{2 \Delta x} = 0 :label: adv_1d_dis1 If we rearrange equation :eq:`adv_1d_dis1` we finally have: .. math:: u_{i}^{n+1}=u_{i}^{n}+ c \frac{ \Delta t}{2 \Delta x} \Big(u_{i-1}^{n}-u_{i+1}^{n}\Big) :label: adv_1d_dis Stability And Accuracy ^^^^^^^^^^^^^^^^^^^^^^ Let assume that solution of equation :eq:`adv_1d` is periodic in the defined domain. Lets consider :math:`\psi(x,t)` is the solution of equation :eq:`adv_1d`. Since this equation is a linear PDE, solution of this PDE is a sum of Fourier modes: .. math:: \psi(x,t) = \sum_{k=-\infty}^{\infty} a_k e^{ikx} :label: for-mod where :math:`k` is the wave number. Now let us study the stability of solution considering an individual wave number. Let us discretize our space and time domain to :math:`x_{j}=j\Delta x` and :math:`t_{n}=n\Delta t`, for :math:`(j,n=0,1,...)`, then the solution at :math:`x_j` at time :math:`t_n` is :math:`\psi(x_j,t_n)`. We can see how the solution is amplified (or damped) by considering analytical amplification factor :math:`A_a`: .. math:: \psi(x,t_{n+1})= A_a \psi(x,t_{n}) :label: A_a Using this and discretization of the spatial domain and considering the solution of one wave number, we have .. math:: \psi_{j}^{n}=\psi(x_j,t_n) = A_a^n e^{ik x_j} = A_a^n e^{i \Delta x} :label: wave Using the same idea for numerical scheme, we can also consider numerical amplification factor :math:`A`: .. math:: \Psi_{j}^{n}=\Psi(x_j,t_n) = A^n e^{ik x_j} = A^n e^{i \Delta x} :label: wave-n where here :math:`\Psi` is our aproximated solution using a linear numerical scheme. Now we can use this equation in our numerical scheme and analyse its stability. Using equation :eq:`wave` in equation :eq:`adv_1d_dis` we have .. math:: A= 1 - i \frac{ \Delta t}{ \Delta x} sin(k \Delta x) :label: FTCS_A In order that our numerical method to be stable, :math:`|A|` should be less than 1. Calculating :math:`|A|` we have .. math:: |A|^2 = 1 + ( \frac{ \Delta t}{ \Delta x} sin(k \Delta x) )^2 :label: FTCS_A2 The above equation can never be able to be less that or equal to 1. Practically this method is not usefull since using von Neumann Stability Analysis, FTCS method is unconditionally unstable. Upwind and Downwind Methods --------------------------- Upwind and downwind methods refer to those methods that the spatial differences are skewed in the flow direction. The simplest upwind and downwind methods are the discribed by backward (:math:`c > 0`) or forward (:math:`c < 0`) spatial difference and the temporal forward Euler methods, respectively. As it is mentioned, the choice between the two spatial methods are dictated by the sign of :math:`c`, meaning the flow direction. Using the descritization methods described above we get: .. math:: \frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t} + \gamma_{min} \Big(\frac{u_{i+1}^{n}-u_{i}^{n}}{ \Delta x} \Big)+ \gamma_{max} \Big(\frac{u_{i}^{n}-u_{i-1}^{n}}{ \Delta x} \Big) = 0 :label: adv_1d_upwind1 where :math:`\gamma_{min}:=min(a,0)` and :math:`\gamma_{max}:=max(a,0)`. Rearranging equation :eq:`adv_1d_upwind1` we have: .. math:: u_{i}^{n+1}= \alpha_{max} u_{i-1}^{n}+(1+\alpha_{min}-\alpha_{max})u_{i}^{n} -\alpha_{min}u_{i+1}^{n} :label: adv_1d_upwind where :math:`\alpha_{min}:=\frac{\gamma_{min} \Delta t}{\Delta x}` and :math:`\alpha_{max}:=\frac{\gamma_{max} \Delta t}{\Delta x}`. Considering the solution is periodic, we can rewrite equation :eq:`adv_1d_upwind` in its matrix form: .. math:: \mathbf{U}^{n+1}= \mathbf{A}^{n+1}\mathbf{U}^{n} :label: mat_eq where :math:`\mathbf{U}^{n+1}=[u_{0}^{n+1},u_{1}^{n+1},...,u_{N-1}^{n+1},u_{N}^{n+1}]^{T}` and :math:`\mathbf{U}^{n}=[u_{0}^{n},u_{1}^{n},...,u_{N-1}^{n},u_{N}^{n}]^{T}` and matrix :math:`\mathbf{A}` for equation :eq:`adv_1d_upwind` is: .. math:: :label: matrixA \small\mathbf{A}=\small\begin{bmatrix} (1+\alpha_{min}-\alpha_{max}) & -\alpha_{min} & 0 & 0 & 0& ... & \alpha_{max} \\ \alpha_{max} & (1+\alpha_{min}-\alpha_{max}) & -\alpha_{min}& 0 & 0 &... & 0 \\ 0 & \alpha_{max} & (1+\alpha_{min}-\alpha_{max}) & -\alpha_{min}& 0 &... & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots & \vdots \\ 0 &0 &... & \alpha_{max} & (1+\alpha_{min}-\alpha_{max}) & -\alpha_{min} & 0 \\ 0 &0 &... & 0 & \alpha_{max} & (1+\alpha_{min}-\alpha_{max}) & -\alpha_{min} \\ -\alpha_{min} &0 &... & 0 & 0& \alpha_{max} & (1+\alpha_{min}-\alpha_{max}) \\ \end{bmatrix}_{N\times N } As we mentioned above, if (:math:`c > 0`) we refer the discretization methods above as upwind and if (:math:`c < 0`) as downwind. Stability And Accuracy ^^^^^^^^^^^^^^^^^^^^^^ Using the same idea from previous section, we calculate the :math:`A` for upwind method: .. math:: A = (1-C)+C e^{-i k \Delta x} :label: Upwind_AUp where :math:`C=\frac{c \Delta t}{\Delta x}`. Now calculating its norm, we have: .. math:: |A|^2 = 1 - 2 C (1-C)(1-cos(k \Delta x)) :label: Upwind_A If this method is stable, then :math:`|A| \leq 1`. Then we have, .. math:: 1 - 2 C (1-C)(1-cos(k \Delta x)) \leq 1 \to 0 \le C (1-C)(sin^2 (k \Delta x / 2)) :label: Upwind_A2 which leads to .. math:: C= \frac{c \Delta t}{\Delta x} \leq 1 :label: cfl_cond Thus this method is conditionally stable. Using the same analogy, for the downwind method we have .. math:: |A|^2 = 1 - 2 C (1+C)(1-cos(k \Delta x)) :label: dwind_A and if the method is stable, then :math:`|A| \leq 1` leading to .. math:: C= \frac{c \Delta t}{\Delta x} \geq -1 :label: cfl_cond-dw Python Code (Matrix form) ^^^^^^^^^^^^^^^^^^^^^^^^^ Here is a python code for modeling the 1D linear advection equation using upwind method described above. .. literalinclude:: python_codes/LA1D.py :language: python .. image:: images/LA1D_Python.png :width: 500px :height: 400px :alt: alternate text :align: center For the explanaition of the code, please take a look at the youtube video. Matlab Code (Matrix form) ^^^^^^^^^^^^^^^^^^^^^^^^^ Here is a Matlab code for modeling the 1D linear advection equation using upwind method described above. First we define a class called LinearAdvection1D: .. literalinclude:: matlab_codes/LinearAdvection1D_upwind.m :language: matlab For the main code we use LinearAdvection1D class in the folloing m file. .. literalinclude:: matlab_codes/LA1D_upwind.m :language: matlab For the explanaition of the code, please take a look at the youtube video. .. image:: images/LA1D_Matlab.png :width: 500px :height: 400px :alt: alternate text :align: center Julia Code (Matrix form) ^^^^^^^^^^^^^^^^^^^^^^^^^ Here is a Julia code for modeling the 1D linear advection equation using upwind method described above. .. literalinclude:: julia_codes/LinearAdvection1D-upwind.jl :language: julia For the explanaition of the code, please take a look at the youtube video. .. image:: images/LA1D_Julia.png :width: 500px :height: 400px :alt: alternate text :align: center Modified Equation ^^^^^^^^^^^^^^^^^ Modified equation first was introduced by Warming et al. which is a powerful tool of analysis the truncation error of time-dependent linear PDEs. In particular, the effect of diffusion or dispersion of the error terms can be studied using modified equation. Considering equation :eq:`adv_1d_upwind`, the modified equation refers to the following formal equation: .. math:: v(x,t+\Delta t) = \alpha_{max} v(x-\Delta x,t) + (1+\alpha_{min}-\alpha_{max}) v(x,t) - \alpha_{min} v(x+\Delta x,t) :label: modif-upwind Assuming :math:`v` is :math:`C^{\infty}`, one can use taylor expansion of all terms in the above equation and deduce the folloing: .. math:: v_{t} +\Delta t v_{tt} +O(\Delta t^2)= -(\gamma_{max}+\gamma_{min}) v_x + \frac{\Delta x}{2}(\gamma_{max}-\gamma_{min}) v_{xx} - \frac{\Delta^{2} x}{6}(\gamma_{max}+\gamma_{min}) v_{xxx} +O(\Delta x^3) :label: modif-upwind1 using :math:`\gamma_{max}+\gamma_{min}= c` and :math:`\gamma_{max}-\gamma_{min}= sign(c)c=|c|` we have: .. math:: v_{t}+ c v_x =-\Delta t v_{tt} + \Big(\frac{\Delta x}{2} |c|\Big)v_{xx} - \Big(\frac{\Delta^{2} x}{6}c\Big) v_{xxx} +O(\Delta^2 t,\Delta^3 x) :label: modif-upwind2 taking the partial derivative of equation :eq:`modif-upwind2` with respect :math:`t` and :math:`x` we have: .. math:: v_{tt}+ c v_{xt} =-\Delta t v_{ttt} + \Big(\frac{\Delta x}{2} |c|\Big)v_{xxt} - \Big(\frac{\Delta^{2} x}{6}c\Big) v_{xxxt} +O(\Delta^2 t,\Delta^3 x) :label: modif-upwind3 .. math:: v_{tx}+ c v_{xx} =-\Delta t v_{ttx} + \Big(\frac{\Delta x}{2} |c|\Big)v_{xxx} - \Big(\frac{\Delta^{2} x}{6}c\Big) v_{xxxx} +O(\Delta^2 t,\Delta^3 x) :label: modif-upwind4 If we multiply equation :eq:`modif-upwind4` by :math:`-c` and sum it with equation :eq:`modif-upwind3` we get: .. math:: v_{tt} =c^2 v_{xx}+\Delta t\Big(c v_{ttx} -t v_{ttt}\Big) + \Delta x \Big(\frac{|c|}{2} v_{xxt} - c \frac{|c|}{2} v_{xxx} \Big) - \Delta^2 x \Big(\frac{c }{6} v_{xxxt} +\frac{c^{2} }{6} v_{xxxx}\Big) +O(\Delta^2 t,\Delta^3 x) :label: modif-upwind5 Substituting equation :eq:`modif-upwind5` to equation :eq:`modif-upwind2` we have: .. math:: v_{t}+ c v_x = \frac{c \Delta x}{2} \Big(Sign(c) - \frac{c \Delta t}{\Delta x} \Big)v_{xx} +O(\Delta^2 t,\Delta^3 x) :label: modif-upwind-f Formally, :eq:`modif-upwind-f` is the PDE which actually solved by upwind method. As can be seen, the first term in the right hand side creates artificial diffusion for equation :eq:`adv_1d` when solved numerically using upwind method. This shows that the upwind method creates numerical diffusion as wave propagates during time. This equation also shows that the method is first order accurate. As we mentioned in equation :eq:`cfl_cond`, upwind method is stable if :math:`|\frac{c \Delta t}{\Delta x}| \le 1`. This means that the diffusion coefficient in equation :eq:`modif-upwind-f` :math:`\frac{c \Delta x}{2} \Big(Sign(c) - \frac{c \Delta t}{\Delta x}\Big)` is positive. Error due to Diffusion and Dispersion ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Similar to von-Neumann stablity analysis, assume a wave-like solution of .. math:: \psi(x,t) = e^{i(kx-\omega t)} = e^{-i \phi_a} :label: wave-ds where :math:`\phi_a` is the wave phase. Substituting equation :eq:`wave-ds` to :eq:`adv_1d`, we get :math:`\omega (k)=ck` which corresponds to the relation between each wave frequency and its wavenumber :math:`k`. In general this relation is called dispersion. if we use discretization of :math:`x_{j}=j\Delta x` and :math:`t_{n}=n\Delta t`, for :math:`(j,n=0,1,...)`, then the solution at :math:`x_j` at time :math:`t_n` is :math:`\psi(x_j,t_n)`. We can now consider the analytical amplification factor: .. math:: A_a =\frac{\psi(x_j,t_{n+1})}{\psi(x_j,t_{n})}=e^{-i\omega \Delta t} :label: A_a-n Then wave phase also can be calculate using :math:`\omega (k)=ck` we have: .. math:: \phi_a = \omega \Delta t = c k \Delta t= CFL \times k \Delta x :label: phi_a Using equations :eq:`A_a-n` and :eq:`wave-ds`, one can show that: .. math:: \psi(x_j,t_n) = \psi_0 e^{i(kx-\omega t)} = \psi_0 A_a^n e^{i k x_j} :label: wave_disc Since :math:`\omega (k)=ck` for linear advection equation, we can calculate the phase and group velocity, namely: .. math:: V_p =\frac{\omega}{k}=\frac{ck}{k}=c :label: phase_v_ex .. math:: V_g =\frac{\partial \omega}{\partial k}=c :label: group_v_ex This shows that the phase and group velocity are the same, meaning that the wave is non-dispersive. Now consider the upwind method and using the same termiology of equations :eq:`wave_disc` and :eq:`A_a-n`, we have: .. math:: \Psi(x_j,t_n) = \Psi_0 A^n e^{i k x_j} :label: wave_disc_upw where :math:`A` is the numerical amplification factor and is defined as: .. math:: A = \frac{\Psi_j^{n+1}}{\Psi_j^{n}} :label: A Note that :math:`A_a` is derived from the analytical solution of the PDE and that :math:`A` is derived from numerical approximation of the same PDE, which might not be in perfect agreement. As we seen before, in the von-Neumann analysis of upwind method, :math:`A` was calculated in equation :eq:`Upwind_AUp`. We can compare :math:`A_a` and :math:`A` by considering the ratio their norm: .. math:: \epsilon_A = \frac{|A|}{|A_a|}= \sqrt{1-4C(1-C)sin^2(\frac{k \Delta x}{2} )}=\sqrt{1-4C(1-C)sin^2(\frac{\chi}{2} )} :label: difer .. image:: images/ED.png :width: 500px :height: 400px :alt: alternate text :align: center In the figure above we can see for different values of :math:`C`, how :math:`\epsilon_A` changed as the frequency increases. As can be seen for :math:`C=1`, both analytical and numerical amplification magnitude are the same and equal to 1. This is the case that numerical solution is the exact solution of the advection equation. However, as :math:`C` decreases, for higher frequencies, :math:`\epsilon_A` decreases dramatically, meaning that the numerical amplitude is smaller than the exact one. We can also do the same analysis for the phase which is associated to dispersion. Defining :math:`\epsilon_\phi` as the ratio between the numerical and analytical phase, and using equation :eq:`Upwind_AUp`, we have .. math:: \epsilon_\phi = \frac{|\phi|}{|\phi_a|}= \frac{1}{C \chi} \frac{C sin(\chi)}{1-C (1-cos(\chi)} :label: phier .. image:: images/EPHI.png :width: 500px :height: 400px :alt: alternate text :align: center In the figure above we can see for different values of :math:`C`, how :math:`\epsilon_\phi` changed as the frequency increases. As can be seen, there is no dispersion error for :math:`C=1` and :math:`C=0.5`. for the values of :math:`C<0.5`, :math:`\epsilon_\phi<1` meaning that the numerical solution moves slower than the physical one. On the other hand for :math:`C>0.5`, this is vice versa. Method of lines ^^^^^^^^^^^^^^^^ The numerical method of lines (MOL) is a technique for solving PDEs. In this method we discretize all but one dimension (for example temporal terms) and then integrate the semi-discrete problem as a system of ODEs or DAEs. For using this method, the PDE problem should be well posed as an initial value problem in at least one dimension, since the ODE and DAE integrators used are initial value in problem solvers. For this reason MOL cannot be used directly on purely elliptic partial differential equations, such as Laplace's equation. However for solving Laplace's equation using MOL, "method of false transients" can be applied or "semi-analytical method of lines" can be used. For simplicity, let us assume equation :eq:`adv_1d` with :math:`c>0`. Then the semi discretization of this equation using backward spatial descritization is: .. math:: u_t =- c \frac{u_{i}^{n}-u_{i-1}^{n}}{ \Delta x} :label: semi-d We have wide range of ODE solvers in Matlab, Julia and python, which we use some of them here. Python Code (Using MOL) ^^^^^^^^^^^^^^^^^^^^^^^^^ Here is a python code for modeling the 1D linear advection equation using MOL. .. literalinclude:: python_codes/Semi_LA1D_Python.py :language: python For the explanaition of the code, please take a look at the youtube video. .. image:: images/semi-LA1D_Python.png :width: 500px :height: 400px :alt: alternate text :align: center Matlab Code (Using MOL) ^^^^^^^^^^^^^^^^^^^^^^^^^ Here is a Matlab code for modeling the 1D linear advection equation using MOL. .. literalinclude:: matlab_codes/Semi_LA1D_Matlab.m :language: matlab For the explanaition of the code, please take a look at the youtube video. .. image:: images/semi-LA1D_Matlab.png :width: 500px :height: 400px :alt: alternate text :align: center Julia Code (Using MOL) ^^^^^^^^^^^^^^^^^^^^^^^^^ Here is a Julia code for modeling the 1D linear advection equation using MOL. .. literalinclude:: Julia_codes/Semi_LA1D_Julia.jl :language: julia For the explanaition of the code, please take a look at the youtube video. .. image:: images/semi-LA1D_Julia.png :width: 500px :height: 400px :alt: alternate text :align: center Lax-Wendroff Method ----------------------- Lax-Wendroff method belongs to the class of conservative schemes. As we can see later, this method provide a second order (in time and space) numerical solution. There are many ways to derive this method. We use one of these methods. Our main aim is to find the numerical solution at :math:`t=t_n+1` knowing solution at :math:`t=t_n`. For this reason we write the taylor expansion of :math:`u(x,t+\Delta t)` with respect to time, namely: .. math:: u(x,t+\Delta t) = u(x,t) + \Delta t u_{t}(x,t) + \frac{\Delta^2 t}{2} u_{tt}(x,t) + O(\Delta^3 t) :label: taylor-lax But from equation :eq:`adv_1d` we know that :math:`\partial_{t}=-c\partial_{x}`, thus we can substitude this expression in equation :eq:`taylor-lax` and we get: .. math:: u(x,t+\Delta t) = u(x,t) -a \Delta t u_{x}(x,t) +a^2 \frac{\Delta^2 t}{2} u_{xx}(x,t) + O(\Delta^3 t) :label: taylor-lax1 Now we can spatially approximate equation :eq:`taylor-lax1` using central difference scheme for :math:`u_{x}(x,t)` and :math:`u_{xx}(x,t)` which results in: .. math:: u_{i}^{n+1} = a_{-1} u_{i-1}^{n} + a_{0} u_{i}^{n} +a_{1} u_{i+1}^{n} :label: lax-w where :math:`a_{-1}=\frac{1}{2}C(1+C)`, :math:`a_{0}=1-C^2`, :math:`a_{1}=\frac{-1}{2}C(1-C)` and :math:`C=\frac{c \Delta t}{\Delta x}`. Similar to upwind method, we can rewrite equation :eq:`lax-w` in its matrix form .. math:: :label: matrixA-lax \small\mathbf{A}=\small\begin{bmatrix} a_{0} & a_{1} & 0 & 0 & 0& ... & a_{-1} \\ a_{-1} & a_{0} & a_{1}& 0 & 0 &... & 0 \\ 0 & a_{-1} & a_{0} & a_{1}& 0 &... & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots & \vdots \\ 0 &0 &... & a_{-1} & a_{0} & a_{1} & 0 \\ 0 &0 &... & 0 & a_{-1} & a_{0} & a_{1} \\ a_{1} &0 &... & 0 & 0& a_{-1} & a_{0} \\ \end{bmatrix}_{N\times N } Stability And Accuracy ^^^^^^^^^^^^^^^^^^^^^^ Using the same idea from FTCS section, we calculate the :math:`A` for Lax-Wendroff method: .. math:: A = a_{0}+ a_{-1} e^{-ik \Delta x}+ a_{1} e^{ik \Delta x} :label: Lax_AU Now calculating its norm, one can show that we have: .. math:: |A|^2 = 1 - 4 C^2 (1-C^2)sin^4(\frac{\chi}{2}) :label: Lax_A2 Where :math:`\chi=k \Delta x`. The method is stable if :math:`|A| \leq 1` which leads to .. math:: |\frac{c \Delta t}{\Delta x}| \leq 1 :label: cfl_cond-lax Thus this method is conditionally stable. Python Code (Using Matrix form) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ For the code, we just need to modify the ::UpwindMatrixAssembly:: the 1DLinearAdvection class and also in the main code in the upwind Python code. .. literalinclude:: python_codes/LA1D_LaxWendroff-Matrix.py :language: python .. literalinclude:: python_codes/LA1D_LaxWendroff-Solve.py :language: python For the explanaition of the code, please take a look at the youtube video. .. image:: images/LA1D_laxWendroff_Python.png :width: 500px :height: 400px :alt: alternate text :align: center Matlab Code (Using Matrix form) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ For the code, we just need to modify the ::UpwindMatrixAssembly:: the 1DLinearAdvection class and also in the main code in the upwind Matlab code. .. literalinclude:: matlab_codes/LA1D_LaxWendroff-Solve.m :language: matlab .. literalinclude:: matlab_codes/LA1D_LaxWendroff-Matrix.m :language: matlab For the explanaition of the code, please take a look at the youtube video. .. image:: images/LA1D_laxWendroff_Matlab.png :width: 500px :height: 400px :alt: alternate text :align: center Julia Code (Using Matrix form) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ For the code, we just need to modify the ::UpwindMatrixAssembly:: the 1DLinearAdvection type and also in the main code in the upwind Julia code. .. literalinclude:: julia_codes/LA1D_LaxWendroff-Matrix.jl :language: julia For the explanaition of the code, please take a look at the youtube video. .. image:: images/LA1D_laxWendroff_Julia.png :width: 500px :height: 400px :alt: alternate text :align: center Modified Equation ^^^^^^^^^^^^^^^^^ Similar to previous section, considering equation :eq:`lax-w`, the modified equation refers to the following formal equation: .. math:: v(x,t+\Delta t) = a_{-1} v(x-\Delta x,t) + a_{0} v(x,t) + a_{1}v(x+\Delta x,t) :label: lax-w-modified1 Similar to previous section one can write the taylor expansion of each term and replace the temporal term by spatial terms and end up with the following: .. math:: v_{t}= \frac{\Delta x}{6}(C^2 - 1) v_{xxx}+ O(\Delta^3 x,\Delta^3 t) :label: lax-w-modified This shows that lax-Wendroff is second order accurate in time and space. Note that there is no diffusion term and only dispersion term can be seen here. Error due to Diffusion and Dispersion ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Similar to previous section we can analysis the behavior of numerical solution by analysing :math:`\epsilon_A` and :math:`\epsilon_\phi`. As we seen before, in the von-Neumann analysis of upwind method, :math:`A` was calculated in equation :eq:`Lax_AU`. Thus we can calculate :math:`\epsilon_A`: .. math:: \epsilon_A = \frac{|A|}{|A_a|}= \sqrt{1 - 4 C^2 (1-C^2)sin^4(\frac{\chi}{2})} :label: difer2 .. image:: images/ED-Lax-Wendroff.png :width: 500px :height: 400px :alt: alternate text :align: center As can be seen here,