Using Anaconda ipython

 Start Anaconda command prompt and type in the following command :
                                                  ipython
This starts the ipython console which is a pretty neat terminal. Set of imports can be done using a single command in the ipython console:
                                               %pylab
But to make life easier for array operations, incorporating mathematical functions it is always best to start ipython using the command in the command prompt:
                                            ipython --pylab

For more information regarding the above check :
https://stackoverflow.com/questions/20525459/ipython-pylab-vs-ipython


   
             

Problems faced in Spyder Python

Problems faced in Python Spyder

1. Spyder faced problem of failure to start Kernel:
   This was solved by using the following command :
                 conda install -c anaconda python  
   After this command spyder was started from the Anaconda terminal using command:
                               spyder
   Now the problem was finally resolved using resetting Spyder to defaults using:
                          tools -> Reset Spyder to factory defaults
2.  There was a problem that was accompanied using the above method, which is the installation of the latest version of PyQt5 which overrides PyQt4. Hence to force the downgrade to PyQt4 the following command was used:         
   Ref :https://github.com/jnicoleoliveira/SPECData/issues/6

                 conda install --channel https://conda.anaconda.org/conda-forge pyqt
                 conda install -c anaconda pyqt=4.11.4

PCA - Post 4

This post is going to extend on the topic scaling introduced to you in the last post.
 I had given you a statement that scaling using error variances actually is the best scaling technique possible , although I hadn't given any proofs , the examples that has been demonstrated in the assignments stands for itself .

After all the proof of the pudding is in the eating . One more advantage this scaling gives is that we can precisely find the rank of the data matrix , using this type of a scaling . An intuition regarding this can be given as follows , since the last few eigenvalues of the data matrix are are supposed to be zero , it is only the errors that make it non zero therefore scaling the data matrix using error variances should make the last set of eigenvalues = 1 . This is simply a hand waving or may be a even worse kind of an explanation , but it can be proved ,that the above statement is true .

Now that I had made some statements , I present to you the algorithm for the MLPCA briefly introduced to you in the last post.Although how we find the error variances , is yet to be discussed , I had given you one particular method for finding the error variances but such a method simply puts more load on the experimenter side which most people wouldn't actually like , so is there any method for actually finding the error variances from the data itself? It turns out that yes there is a method to do it , it is called the IPCA. 

MLPCA

Algorithm : Let \( \Sigma  _{e}\) be the error covariance matrix of the measurements , $ L^{'}L = \Sigma  _{e} $ . If the error covariances turn out to be zero then $\Sigma_{e}$ is nothing but a diagonal matrix , the following algorithm will do nothing but scaling .
Now using $L$ , we get
$Z_{s} = L^{-1} Z$, Now applying SVD to $Z_{s}$ , we get
$Z_{s} = U_{s} S_{s} V_{s} ^{'}$

Now taking the last few vectors of $C_{s} = U_{s}^{( : ,(k+1) : end)}$  corresponding to singular values = 1 , we actually get the constraint set of equations with respect to the scaled matrix $Z_{s}$ . Where k is the row rank of the data matrix.

Now to converting the constraint set to original matrix we get  :

$C = L^{-1}C_{s}$ . Here C is the constraint set corresponding to the original matrix .

Now , still we hadn't proposed the IPCA , which will actually help you in finding the error covariance matrix with a catch , which we shall discuss in the following blog posts.

The Auto Scaling algorithm 

It is similar to the MLPCA algorithm but instead of  $\Sigma  _{e}$ being the error covariance matrix , it is simply the covariance matrix of the data matrix .

PCA - Post 3

The previous two posts would have given you an intuitive understanding of the mathematics involved in PCA . This post is going to present to you dimensionality reduction aka Data Compression.
Before we learn about Data compression ,we need to discuss few terminologies

Few Terminologies

1. Scores Matrix : SVD of Z gives you $U , S , V$  where $U$ is called the scores matrix.
                               The scores matrix need not be the complete set of $U$ ,it can be the set of vectors corresponding to the first largest eigenvalues.
2. Loadings Matrix : Similarly $V$ matrix is called the loadings matrix.

Dimensionality Reduction

With the above terminologies we will now use them to obtain lower dimensional data . 

We had already proved that in the previous post that $U$ matrix consists of set of vectors of maximum variance . Now projecting the data (mathematically known as Dot product ) onto these directions give the lower dimensional data . 

For example let us assume first  5 vectors of the scores matrix  , premultiplying with the data matrix Z gives corresponding scores of the data matrix. Here the $Z_{m*n}$  $m \geq 5 $ , let m = 8 , n=1000 , for this case. Then it can be seen that the size of the data matrix reduces by 3000. 

Matlab Code

scores = U(: ,1:5)' *Z;

This scores can now be used as any other data matrix which can be used for regression .This is called Principal Component Regression (PCR) . 

These outlines the basic concepts of PCA. Now we come to principles dealing with scaling in PCA,this is a highly controversial topic , practices such as auto scaling , or other types of scaling are performed without any proper justification , we shall first discuss auto scaling .   

Scaling in PCA

Auto Scaling : This is nothing but scaling of mean centered data matrix ($Z$) with the standard deviation of the data matrix . The performance of this type of scaling is usually better than no scaling ,but it mostly depends on the error variances of the data matrix. 

It can be proved using first principles that scaling using error variances of the data matrix performs best . Since obviously the variances of the data matrix need not be proportional to the data matrix at all , there is no guarantee that auto scaling needs to work well . It is to be seen that it performs very bad at times .

MLPCA :

Now that I have claimed that scaling using the error variances , the questions is how to actually find the error variances? 

This can be found using repeated experiments of the same parameters . For example let us talk about an experiment in which we have to measure the absorption spectra for a wide variety of wavelengths for a particular solution (containing various elements ) . Then if for the same a solution of fixed concentration we perform repeated experiments we can actually find the error variances , hence scale the data matrix . Performing PCA after this type of scaling is called MLPCA. 

The mathematical details , algorithms will be discussed in the following posts since it is a vast topic in itself.

PCA - Post 2

 Although I had mentioned about the number of constraint equations based on the row rank of the matrix , when it comes to real life scenario we don't know the actual rank of the matrix since the data is always subject to errors, hence although the data is obtained from an experiment constrained by 3 equations , the row rank of the data matrix will most probably be equal to the number of rows due to errors in measurements , thereby mathematically we can say that it isn't constrained by nay equations when it is actually not the case.

 But heuristically we can actually find the number of constraints approximately by many methods .

First we will look at a different approach to derive the same constraint matrix using eigen value decomposition .

A Different Approach

$Z = Z_{m*n}$

$ZZ^{'} = (USV^{'} )*(VS^{'}U^{'} )$

$\implies ZZ^{'} = (USS^{'}U^{'} )$

Hence , to obtain the constraint matrix , only the U matrix is necessary .

So the new Matlab code : 

[U, Sdash] = eig($ZZ^{'}$)  , where Sdash = $SS^{'}$

Here U, Sdash are the eigenvectors , eigenvalues respectively of the matrix $ZZ^{'}$  . 

Hence it can clearly be seen that the singular values of the data matrix A is equal to the square root of the eigenvalue of the matrix $ZZ^{'}$ .  

The eigenvalues actually give a sense of the statistical variance of the equations whose coefficients are given by the U (eigenvectors) .  For example let $U_{1}$ be the coefficients of an equation , then the statistical variance of that equation  is given by $U_{1}^{'} \Sigma U_{1} $ ,where $\Sigma $ is the covariance matrix across the variables/features  .

Before we actually prove that the eigenvalues of the $ZZ^{'}$ (matrix same as covariance matrix differing by a scale)  actually give the variance , we shall for now accept it as a fact and move forward to predict the number of constraints using it .We can now by using the eigenvalues approximately identify the number of constraints ,how do we do that ?
Let us take an example where in an equation whose coefficients are given by $U_m$ , have an eigenvalue = $10^{-3}$ ,then what this means is that the statistical variance of that equation  $\approx 0$ .
Using the same principle, some methods (although very loose methods) are introduced ,rigorous methods will be introduced in the following posts.

Preliminary methods 

1. SCREE Plot
   This basically is a plot of the eigen values of the covariance matrix (which can differ by a scale)  . This plot basically gives you a visual knee , which can be then used as a heuristic to find the number of constraints . Hence heuristically we can see that the number of constraints is approximately equal to 7. 



Image Courtesy : http://ba-finance-2013.blogspot.in/2012/09/scree-plots-interpretation-and.html



2. Variance captured 
    At times people tend to also use the amount of variance captured  , to find the number of constraints .  the formula for the variance captured is given by:
    $ VarCapt = \frac{\sum _ {i =1} ^{u} e_{i}}{\sum _ {i =1} ^{m} e_{i}}$,where $ u \leq m $
    Here m is the same as the row size of the data matrix introduced in the previous blog. 
Here the denominator is the sum of all eigenvalues , the numerator is the sum of the u largest eigenvalues .
Hence keeping some threshold say 95% for the variance captured (VarCapt) we can find the value of u , which is nothing but a measure of the row rank of the data matrix.  

In the above methods I had briefly introduced some very basic methods to approximately find the number of constraints , without giving proof of how eigenvalues of the covariance matrix actually gives the sense of the statistical variance of the equation whose coefficients are given by the corresponding eigenvector . So let us actually look at the proof :

Proof 

Another way of looking at PCA , is to find a set of unit orthogonal vectors acting as coefficients of an equation to give maximum variances . Hence by posing this as an optimization problem we get 

max $(w^{'} \Sigma w) $ ,such that $w'w = 1$ ,where $w$ is the coefficient vector of the equation , $\Sigma$ the covariance matrix.

Now setting up this as a lagrange function we get 

$L = (w^{'}  \Sigma w)  - \lambda(w'w - 1) $ , 

Hence the new optimization problem is max ($L$)

So , now differentiating L we get  

$\frac{\delta L}{\delta w} = 0 = 2 \Sigma w - 2\lambda(w) $

$\implies  \Sigma w = \lambda(w) $  (ie w , $\lambda$ are the corresponding eigenvectors ,eigenvalues of the covariance matrix . 

$\implies  w^{'} \Sigma w = \lambda(w^{'}w) $

But,  

$\frac{\delta L}{\delta \lambda} = 0 = w'w - 1$

$\implies w'w  =  1$

Hence,

$ w^{'} \Sigma w = \lambda $  (proved) 

PCA - An introduction

Introduction 

The origin of PCA is just analogous to TLS (Total Least Squares ) ,used for regression but it is extended to data compression (reducing dimensionality of data), data imputation (imputation of missing data),Iterative PCA (IPCA) (Used for finding the errors in the data )  , Non-Negative Matrix Factorisation  (NMF) (Used for extracting clean data from noisy data - The famous Cock tail party problem), Network Component Analysis (NCA) (Used generally in Bio-Informatics) , finally  Kernel PCA (KPCA) (to model any physical data). PCA is such a powerful tool , which can be used to do all this provided no consideration is taken on the computational power , that is assuming we have enough computational power we can device an algorithm for all the above mentioned .

Introduction to SVD 

PCA uses an mathematical tool in Matlab called a Singular Valued Decomposition (SVD) . This is very much analogous to any other factorisation technique such as the QR , Eigen Value decomposition ,etc . 

Let us define a matrix A , its SVD is given by :

$A_{m*n}  = U_{m*n}S_{n*m}V_{n*m}^{'}$ 

The Properties of the matrix :  The matrix U , V are the left singular ,right singular vectors which are orthonormal to each other , S is a diagonal matrix consisting of singular values (Generic : first few singular values are non zero arranged in descending order, followed by zero singular values). 

(ie $s_{1} \geq s_{2} \geq s_{3} \cdots  s_{r} \geq s_{r+1} \approx 0$ .(Here singular values after the rth value are approximately equal to 0).

It can now be easily seen that :

$AA^{'} = (USV^{'} )*(VS^{'}U^{'} ) $

$\implies AA^{'} = U(SS^{'})U^{'}$

Hence it can be seen that $ U $, $SS^{'}$ is set of eigen vectors ,eigen values respectively of the matrix $AA^{'}$  . It can also be noted that if  $A$ is a data matrix that contains each data points as a column vector (ie the matrix sizing is Variable x No of datapoints ) , then the matrix  $AA^{'}$ kind of acts like a covariance matrix across the variables . 

Similarly ,

$A^{'}A = (VS^{'}U^{'} )*(USV^{'} ) $

$\implies A^{'}A = V(S^{'}S)V^{'}$

Hence it can be seen that $ V $, $S^{'}S$ is set of eigen vectors ,eigen values respectively of the matrix $A^{'}A$  . It can also be noted that if  $A$ is a data matrix that contains each data points as a column vector (ie the matrix sizing is Variable x No of datapoints ) , then the matrix  $AA^{'}$ kind of acts like a covariance matrix across the datapoints. 

Now it can also be seen that 

$A = \sum _{i = 1} ^{k} U_{i}s_{i}V_{i}^{'}$

Where $ U_{i},V_{i}$ is the ith column of the U matrix , V matrix respectively ,$s_{i}$ represents the ith non zero singular value  ,where $ k = min(m,n) $ 

Now that adequate details have been given about the SVD , the properties of the decomposed matrices  , I hope you can yourself formulate the method that is required to perform the decomposition (ie SVD) .

Algorithm 

A setup of the problem :

For all purposes let us define some basic conventions that will be used through out the series of blogs under MVDA , 

Let A  - $A_{m*n}$be the data matrix (Variable x No of Datapoints ) .

Then matrix $Z = A - \bar{A}$ , here $\bar{A}$ is the mean of the data matrix .

for all purposes we shall now be dealing with Z matrix. 

Now some inferences we shall make on Z matrix :

First let us assume the Z matrix to be a matrix with no errors , the system following an exact linear relation then the number of constraints  = m - row rank of matrix (Z).  

Algorithm Approach  :

Then , let us perform an SVD on Z ,

then if we take the eigen vector corresponding ($U_{m}$) to the last singular value gives we get  :
Note : when $m \leq n$ (a well posed data matrix) then $k = m$ ,hence : 

$ U_{k}^{'}  Z = U_{k}^{'}* (\sum _{i = 1} ^{k} U_{i}s_{i}V_{i}^{'})$

$\implies U_{k}^{'}  Z =  s_{k}V_{k}^{'} \approx  0$. (pre-multiply with the vector)

Now you may see how the TLS algorithm discussed in my previous blog even originated .

If the data matrix is ill posed $m \geq n$ then $m \geq k$ hence

$ U_{m}^{'}  Z = U_{m}^{'}* (\sum _{i = 1} ^{k} U_{i}s_{i}V_{i}^{'})$

$\implies U_{m}^{'}  Z =  0$. (pre-multiply with the vector)

Now if there are multiple constraints then we obviously take the last few vectors of the U matrix that can now be called a regression set . 

Now as mentioned that the Z matrix consists of each and every data entry in the form of a column , if the data matrix that you are using is the other way around the you can post multiply with the last few vectors of $V$ corresponding to least singular values .

Matlab Function

Since we have seen a brief intro to SVD  , I shall now introduce to you the matlab function that would do the decomposition for you :

Matlab Function : [U,S,V] = svd(A).

The above function can always be used to do the SVD  , to obtain a regression set.
Since it was already mentioned that $k = min(m,n) $, lets consider a case where $m \leq n$ , then the extra vectors of the $V$ matrix (ie columns of the V matrix from $k+1$ to $n$ columns) are unnecessarily computed , thereby increasing computational power.

Hence we can use a specialized version of the same function under certain conditions:

[U,S,V] = svd(A,'econ')

This version basically ignores the above mentioned columns of the $V$ matrix  , also returns $S$ as a square matrix.
Note : Assuming blog convention of data matrix :
           The econ version is only used when $m \leq n$ .
            If $m \geq n$ then econ shouldn't be used.

Lid Driven Cavity - Post 3

Code Files$^1$

1. Main.m (The file to be run)
2. streamfunc.m (Solves the stream vorticity function)
3. velocity.m (Solves the velocity function)
4. omega.m (Iterates the vorticity transport equation)
5. BC.m (Implements the Boundary conditions)
6. Upwind.m (Implements the upwind scheme)
7. rmse psi.m (Computes Errors)

Terminologies used in code$^2$

1. w = $\omega$(Vorticity)
2. psi = $\psi$(StreamFunction)
3. u = (U −V elocity)
4. v = (V −V elocity)
5. u0 = $u_0$(Lid Velocity)
6. psi_1 = $c_1$
7. gamma = $\gamma$ (Kinematic viscosity)
8. dx = ∆x
9. dy = ∆y
10. alpha = Relaxation parameter for stream vorticity function
11. alpha1 = Relaxation parameter for vorticity transport equation
12. x = (X −length) of cavity
13. y = (Y −length) of cavity
14. Nx = $N_x$(No of grid points along x-direction
15. Ny = $N_y$(No of grid points along y-direction
16. err - Percentage error in $\psi$
17. err1 - Percentage error in $\omega$
18. ERR1 - Error in stream vorticity equation
19. ERR2 - Error in vorticity transport equation

-------------------------------------------------------------------------------------------------------------------------
$^1$ Link to code given under section Matlab Code URL
$^2$ Variable Names used in the Matlab code implementing the Stream Vorticity approach for the Modified Lid Driven Cavity Problem

Lid Driven Cavity - Post 2

(You can find the PDF version of the below blogpost here. )
(Citation : R, AAKASH. "Lid Driven Cavity." Classroom fundaes. Aakash R, 2 May 2017. Web. http://instafundae.blogspot.in/2017/05/lid-driven-cavity-post-2.html.)

Abstract

The lid driven cavity is a standard benchmark problem for any CFD software/Code, Lid driven cavity flows are important in many industrial processing applications such as short-dwell and flexible blade coaters.This problem consists of a 2-D cavity , whose top side moves with  a uniform velocity $u_0$.There are numerous results form the literature for the above problem , which have been used to benchmark the code,its respective results presented in this paper.The benchmarking references used in this paper are A. AbdelMigid [1] ,O. BOTELLA [3],Ghia $et$ $al.$[4]
This paper presents steady state solutions to Modified Lid Driven Cavity(Same as lid driven cavity problem but with inflow and outflow at the left and right hand side boundary) introduced in the book by Hoffman [5]  for low to Middle  $Re$s which can be extended to produce results for the Lid Driven Cavity . It is to be duly noted that the results presented in this paper are just a small part of what the Code (The Matlab code written to solve the Modified LDC) has been tested for,the same code can be used for middle to high level $Re$s too with appropriate grid sizes.
The approach used in solving the above mentioned problem is the standard stream vorticity approach . The numerical schemes used for producing results in this include central ,forward ,backward differencing schemes,Guass seidel PSOR(Point Successive Over Relaxation) methods in solving equation,etc.

The schematic diagram of the Modified Lid Driven Cavity is as follows :

Governing Equations

The governing equations of the stream vorticity approach are Vorticity transport equation,stream function  equation , continuity equations which are coupled .The below equations can be easily derived .
\begin{eqnarray}
\frac{\delta \omega}{\delta t} + u \frac{\delta \omega}{\delta x} +v\frac{\delta \omega}{\delta y} &=& \gamma \left(\frac{\delta ^2 \omega}{\delta x ^2} + \frac{\delta ^2 \omega}{\delta y ^2}\right)\\
  \nabla ^2\psi&=& - \omega\\
u &=&  \frac{\delta \psi}{\delta y}\\
v &=& -\frac{\delta \psi}{\delta x}
\end{eqnarray}
Note the continuity equation is implicitly used in deriving the stream vorticity formulations.

Boundary Conditions
Modified Lid Driven Cavity

Boundary conditions for $u,v,\psi$

Bottom side :
\begin{eqnarray}
u &=& 0\\
v &=& 0\\
\psi &=& c_1
\end{eqnarray}
Top side :
\begin{eqnarray}
u &=& 1\\
v &=& 0\\
\psi &=& c_2
\end{eqnarray}
Right side:
For $1:(j_3 -1)$
\begin{eqnarray}
u &=& 0\\
v &=& 0\\
\psi &=& c_1
\end{eqnarray}
 For $j_3:j_4$
\begin{eqnarray}
\frac{\delta \psi}{\delta x} &=& 0
\end{eqnarray}
 For $(j4+1):end$
\begin{eqnarray}
u &=& 0\\
v &=& 0\\
\psi &=& c_2
\end{eqnarray}

Left side :
 For $1:(j_1 -1)$
\begin{eqnarray}
u &=& 0\\
v &=& 0\\
\psi &=& c_1
\end{eqnarray}
 For $j_1:j_2$
\begin{eqnarray}
\frac{\delta \psi}{\delta x} &=& 0
\end{eqnarray}
 For $(j2+1):end$
\begin{eqnarray}
u &=& 0\\
v &=& 0\\
\psi &=& c_2
\end{eqnarray}

Boundary conditions for $\omega$

\begin{equation}
  \omega  = -\left(  \frac{\delta ^2 \psi}{\delta x^2} + \frac{\delta ^2 \psi}{\delta y^2}\right)
\end{equation}

Discretized Equations

Vorticity Transport Equation

Inner Grid Formulation
Central Differencing Scheme
\begin{equation}
u(i,j)\frac{\omega(i+1,j)-\omega(i-1,j) }{2\Delta x}+v(i,j)\frac{\omega(i,j+1)-\omega(i,j-1) }{2\Delta y} =  \gamma \left( \frac{\omega(i+1,j)-2\omega(i,j)+\omega(i-1,j)}{(\Delta x)^2} + \frac{\omega(i,j+1)-2\omega(i,j)+\omega(i,j-1)}{(\Delta y)^2}\right)
\end{equation}
Rearranging , equating $\Delta x = \Delta y$,we get:
\begin{equation}
\omega(i,j) = \frac{1}{4} \left( \omega(i+1,j)\left(1 - \frac{u(i,j)\Delta x}{2\gamma}\right) + \omega(i-1,j)\left(1 + \frac{u(i,j)\Delta x }{2\gamma}\right) \\ + \omega(i,j+1)\left(1 - \frac{v(i,j)\Delta x }{2\gamma}\right) +  \omega(i,j-1)\left(1 + \frac{v(i,j)\Delta x}{2\gamma}\right)\right)
\end{equation}
Finally applying the PSOR
\begin{equation}
\omega(i,j) =(1-\alpha)\omega(i,j) +  \frac{\alpha}{4} \left( \omega(i+1,j)\left(1 - \frac{u(i,j)\Delta x}{2\gamma}\right) + \omega(i-1,j)\left(1 + \frac{u(i,j)\Delta x }{2\gamma}\right) \\ + \omega(i,j+1)\left(1 - \frac{v(i,j)\Delta x }{2\gamma}\right) +  \omega(i,j-1)\left(1 + \frac{v(i,j)\Delta x}{2\gamma}\right)\right)
\end{equation}
Boundary Formulation

Left Boundary 
\begin{eqnarray}
\omega(1,j )&=& - \left. \frac{\delta ^2 \psi}{\delta x^2} \right \vert _{(1,j)} - \left.\frac{\delta u}{\delta y}\right \vert _{(1,j)}\\
\implies \omega(1,j) &=& 2\frac{\psi(1,j)-\psi(2,j)}{\Delta x^2} - 2\frac{v(1,j)}{\Delta x} -\frac{u(1,j+1)-u(1,j-1)}{(2\Delta y)}
\end{eqnarray}
Right Boundary 
\begin{eqnarray}
\omega(im,j )&=& - \left. \frac{\delta ^2 \psi}{\delta x^2} \right \vert _{(1,j)} - \left.\frac{\delta u}{\delta y}\right \vert _{(1,j)}\\
\implies \omega(1,j) &=& 2\frac{\psi(1,j)-\psi(2,j)}{\Delta x^2} - 2\frac{v(1,j)}{\Delta x} -\frac{u(1,j+1)-u(1,j-1)}{(2\Delta y)}
\end{eqnarray}

Stream Vorticity Equation

Inner Grid Formulation
Central Differencing Scheme
\begin{eqnarray}
\omega = - \left( \frac{\psi(i+1,j)-2\omega(i,j)+\psi(i-1,j)}{(\Delta x)^2} + \frac{\psi(i,j+1)-2\psi(i,j)+\psi(i,j-1)}{(\Delta y)^2}\right)\\
\implies \psi(i,j) = \frac{1}{4}\left(\psi(i,j-1)+\psi(i,j+1)+\psi(i-1,j)+\psi(i+1,j)+ \omega (i,j)\Delta x ^2 \right)
\end{eqnarray}
Finally applying PSOR we get
\begin{equation}
 \psi(i,j) = (1-\alpha)\psi(i,j)+\frac{\alpha}{4}(\psi(i,j-1)+\psi(i,j+1)+\psi(i-1,j)+\psi(i+1,j)+ \omega (i,j)\Delta x ^2 )
\end{equation}

Boundary Formulation

Left Boundary ($j_1 : j_2$):
\begin{eqnarray}
\left. \frac{\delta \psi}{\delta x} \right \vert _{(1,j)} = v
\end{eqnarray}
Using 2nd order accurate scheme we get
\begin{eqnarray}
\psi(1,j) = \frac{4\psi(2,j)-\psi(3,j) + 2\Delta x v(1,j)}{3}
\end{eqnarray}
 Right Boundary ($j_3 : j_4$):
Using 2nd order accurate scheme we get
\begin{eqnarray}
\psi(im,j) = \frac{4\psi(im-1,j)-\psi(im-2,j) - 2\Delta x v(1,j)}{3}
\end{eqnarray}

Velocity Equations

Inner Grid Formulation

Central Differencing Scheme
\begin{eqnarray}
u = \frac{\psi(i,j+1)-\psi(i,j-1)}{(2\Delta  y)}\\
v = \frac{\psi(i+1,j)-\psi(i-1,j)}{(2\Delta  x)}
\end{eqnarray}

Boundary Formulation

Forward Differencing Scheme
\begin{eqnarray}
u &=& \frac{(\psi(j+1,i) - \psi(j,i))}{\Delta y}\\
v &=& 0
\end{eqnarray}

Parameters used for convergence measure
Before we see the algorithm used in solving the problem ,I would introduce few parameters used in checking the convergence of the solution.
The parameters are as follows :

1.  $e_1$: RMSE $\left ( \frac{\psi  ^{k+1} _{i,j} -\psi  ^k _{i,j}}{\psi  ^k _{i,j}} \right )$ (where k is the iteration number)
   Variable name used in the program ("err").
2.  $e_2$ :RMSE$\left ( \frac{\omega  ^{k+1} _{i,j} -\omega  ^k _{i,j}}{\omega  ^k _{i,j}}\right )$ .
   Variable name used in the program ("err1")
3.  $e_3$ : Numerical solution accuracy of Vorticity equation .
   Variable name used in the program ("ERR1")
4. $e_4$ : Numerical solution accuracy of Navier stokes equation .
   Variable name used in the program ("ERR2")

It is to be noted that $e_1<e_2<e_3<e_4$ for the algorithm used in this paper , hence convergence is measured using $e_4$ , divergence measured using $e_1$.

Algorithm

1. Assuming a specific grid size $N_x,N_y $  .
2.  Construct matrices for  $\psi,\omega,u,v$  using the above mentioned size using random values .  
3. Assume appropriate boundary conditions for  $\psi,\omega,u,v$ .
4. Solve for  $\psi$  using the PSOR method using Stream Vorticity function .
5.  Compute  $u,v$  using central difference scheme .
6. Iterate the PSOR method for  $\omega$   using vorticity transport Equation once (using the pre solved  $\psi,u,v $ .
   (Note : Don't solve for  $\omega$ run the PSOR over the entire grid once.)
7. Compute necessary Error Functions, check for convergence.
8. If converged , STOP. Else go to Step 4.

Results

Presented below are few results which can be used for used for benchmarking purposes . 

129x129 Grid 

S.No	U Velocity  (m/s)					V velocity  (m/s)
Distance (m)	Re 100	Re 400	Re 800	Re 1000	Re 3200	Re 100	Re 400	Re 800	Re 1000	Re 3200
1.00	1.0000	1.0000	1.0000	1.0000	1.0000	0.0000	0.0000	0.0000	0.0000	0.0000
0.9760	0.8431	0.7568	0.6799	0.6537	0.5189	-0.0465	-0.0885	-0.1349	-0.1655	-0.3331
0.9687	0.7911	0.6825	0.5955	0.569	0.4747	-0.0621	-0.1224	-0.1885	-0.2311	-0.4441
0.9609	0.7402	0.6155	0.5283	0.5068	0.4609	-0.0777	-0.1577	-0.2434	-0.2967	-0.5218
0.9531	0.6906	0.5567	0.4775	0.4632	0.4583	-0.0930	-0.1936	-0.2974	-0.3586	-0.5537
0.8515	0.2357	0.2883	0.3211	0.3350	0.3465	-0.2388	-0.4452	-0.4241	-0.4061	-0.3679
0.7344	0.0037	0.1603	0.1806	0.1892	0.2000	-0.2134	-0.2522	-0.2295	-0.2356	-0.2300
0.6172	-0.1388	0.0198	0.0483	0.0551	0.0775	-0.0738	-0.0900	-0.0981	-0.0999	-0.1057
0.5000	-0.2087	-0.1159	-0.0705	-0.0620	-0.0367	0.0574	0.0526	0.0295	0.0258	0.0143
0.4531	-0.2134	-0.1722	-0.1161	-0.1089	-0.0800	0.0970	0.1110	0.0798	0.0771	0.0605
0.2813	-0.1570	-0.3244	-0.2982	-0.2808	-0.2405	0.1752	0.2855	0.2758	0.2690	0.2361
0.1719	-0.1013	-0.2376	-0.3564	-0.3852	-0.3434	0.1696	0.2862	0.3515	0.3715	0.3556
0.1016	-0.0641	-0.1418	-0.2456	-0.2923	-0.4297	0.1323	0.2357	0.3070	0.3370	0.4291
0.0703	-0.0464	-0.1001	-0.1783	-0.2236	-0.4036	0.1031	0.1947	0.2625	0.2958	0.4071
0.0625	-0.0418	-0.0897	-0.1610	-0.2029	-0.3798	0.0944	0.1812	0.2475	0.2804	0.3914
0.0547	-0.0371	-0.0792	-0.1435	-0.1817	-0.3499	0.0850	0.1661	0.2302	0.2625	0.3719
0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000

Comparison of my code with Ghia et al (129x129 Grid) 

Grid Points	U-Velocity		Grid Points	V-Velocity
	Present	Ghia et al		Present	Ghia et al
	Re 100			Re 100
129	1	1	129	0	0
126	0.84311	0.84123	125	-0.0621	-0.05906
125	0.79114	0.78871	124	-0.07766	-0.07391
124	0.74015	0.73722	123	-0.09305	-0.08864
123	0.69055	0.68717	122	-0.10817	-0.10313
110	0.23571	0.23151	117	-0.17657	-0.16914
95	0.00375	0.00332	111	-0.23289	-0.22445
80	-0.1388	-0.13641	104	-0.2526	-0.24533
65	-0.2087	-0.20581	65	0.057388	0.05454
59	-0.2134	-0.2109	31	0.178879	0.17527
37	-0.157	-0.15662	30	0.178671	0.17507
23	-0.1013	-0.1015	21	0.164132	0.16077
14	-0.0641	-0.06434	13	0.125833	0.12317
10	-0.0464	-0.04775	11	0.11128	0.1089
9	-0.0418	-0.04192	10	0.103134	0.10091
8	-0.0371	-0.03717	9	0.094379	0.09233
1	0	0	1	0	0

251x251 Grid

S.No	U Velocity  (m/s)				V velocity  (m/s)
Distance (m)	Re 100	Re 400	Re 1000	Re 3200	Re 100	Re 400	Re 1000	Re 3200
1.00	1.0000	1.0000	1.0000	1.0000	0.0000	0.0000	0.0000	0.0000
0.9760	0.8431	0.7568	0.6537	0.5189	-0.0465	-0.0885	-0.1655	-0.3331
0.9687	0.7911	0.6825	0.569	0.4747	-0.0621	-0.1224	-0.2311	-0.4441
0.9609	0.7402	0.6155	0.5068	0.4609	-0.0777	-0.1577	-0.2967	-0.5218
0.9531	0.6906	0.5567	0.4632	0.4583	-0.0930	-0.1936	-0.3586	-0.5537
0.8515	0.2357	0.2883	0.3350	0.3465	-0.2388	-0.4452	-0.4061	-0.3679
0.7344	0.0037	0.1603	0.1892	0.2000	-0.2134	-0.2522	-0.2356	-0.2300
0.6172	-0.1388	0.0198	0.0551	0.0775	-0.0738	-0.0900	-0.0999	-0.1057
0.5000	-0.2087	-0.1159	-0.0620	-0.0367	0.0574	0.0526	0.0258	0.0143
0.4531	-0.2134	-0.1722	-0.1089	-0.0800	0.0970	0.1110	0.0771	0.0605
0.2813	-0.1570	-0.3244	-0.2808	-0.2405	0.1752	0.2855	0.2690	0.2361
0.1719	-0.1013	-0.2376	-0.3852	-0.3434	0.1696	0.2862	0.3715	0.3556
0.1016	-0.0641	-0.1418	-0.2923	-0.4297	0.1323	0.2357	0.3370	0.4291
0.0703	-0.0464	-0.1001	-0.2236	-0.4036	0.1031	0.1947	0.2958	0.4071
0.0625	-0.0418	-0.0897	-0.2029	-0.3798	0.0944	0.1812	0.2804	0.3914
0.0547	-0.0371	-0.0792	-0.1817	-0.3499	0.0850	0.1661	0.2625	0.3719
0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000	0.0000

Comparison with Ghia etal ,and O .BOTELLA:

Grid Points	U-Velocity					Grid Points	V-Velocity
	Present	Ghia et al	Present	Ghia et al	O.BOTELLA		Present	Ghia et al	Present	Ghia et al	O.BOTELLA
	Re 400		Re 1000				Re 400		Re 1000
251	1	1	1	1	1	251	0	0	0	0	0
245	0.75389422	0.75837	0.653659036	0.65928	0.6644227	243	-0.127997501	-0.12146	-0.2311	-0.21388	-0.22792
243	0.67893746	0.68439	0.569419488	0.57492	0.5808359	241	-0.16480199	-0.15663	-0.2967	-0.27669	-0.29369
241	0.61153269	0.61756	0.506788795	0.51117	0.5169277	239	-0.202194664	-0.19254	-0.3586	-0.33714	-0.35532
239	0.5527114	0.55892	0.463245781	0.46604	0.4723329	237	-0.239373072	-0.22847	-0.4134	-0.39188	-0.41038
214	0.29147429	0.29093	0.334964485	0.33304	0.3372212	228	-0.382661909	-0.23827	-0.5226	-0.5155	-0.52644
185	0.16391904	0.16256	0.189150105	0.18719	0.1886747	216	-0.452317228	-0.44993	-0.426	-0.42665	-0.42645
155	0.01930732	0.02135	0.055104972	0.05702	0.0570178	202	-0.383530806	-0.38598	-0.3172	-0.31966	-0.32021
126	-0.1152734	-0.11477	-0.061988317	-0.0608	-0.062056	126	0.052196894	0.05188	0.02584	0.02526	0.0258
114	-0.1730285	-0.17119	-0.108910482	-0.10648	-0.1082	60	0.302110802	0.30174	0.32157	0.32235	0.32536
71	-0.3276167	-0.32726	-0.280799838	-0.27805	-0.28037	58	0.30257526	0.30203	0.33037	0.33075	0.33399
44	-0.2423208	-0.24299	-0.385150406	-0.38289	-0.388569	40	0.282024218	0.28124	0.37335	0.37095	0.37692
26	-0.1428726	-0.14612	-0.292346467	-0.2973	-0.300456	24	0.228366656	0.22965	0.32684	0.32627	0.33304
19	-0.1047669	-0.10338	-0.223618183	-0.2222	-0.222896	21	0.212590596	0.2092	0.30931	0.30353	0.30991
17	-0.0938794	-0.09266	-0.202882065	-0.20196	-0.20233	19	0.200473862	0.19713	0.29582	0.29012	0.29627
15	-0.0829239	-0.08186	-0.181688508	-0.18109	-0.181288	17	0.1868535	0.1836	0.28043	0.27485	0.28071
1	0	0	0	0	0	1	0	0	0	0	0

501x501 Grid:

Note : The Underlined values show some discrepancies with the results of Ghia etal.

Grid Points	U-Velocity			Grid Points	V-Velocity
	Present	Ghia et al	AbdelMigid		Present	Ghia et al	AbdelMigid
	Re 3200	Re 3200	Re 3200		Re 3200	Re 3200	Re 3200
501	1	1	1	501	0	0	0
489	0.518868	0.53236	0.536	485	-0.44412	-0.39017	−0.4106
485	0.474687	0.48296	0.4886	481	-0.52184	-0.47425	−0.5049
481	0.460908	0.46547	0.4661	478	-0.55368	-0.52357	−0.5560
478	0.458339	0.46101	0.461	474	-0.56399	-0.54053	−0.5671
427	0.346469	0.34682	0.3491	454	-0.44061	-0.44307	−0.4462
368	0.200026	0.19791	0.2035	431	-0.37699	-0.37401	−0.3799
310	0.077457	0.07156	0.078	403	-0.31086	-0.31184	−0.3148
251	-0.03674	-0.04272	−0.0369	251	0.014275	0.00999	0.0142
228	-0.07997	-0.86636	−0.0808	118	0.286802	0.28188	0.2864
142	-0.24047	-0.24427	−0.2414	114	0.295331	0.2903	0.2954
87	-0.34339	-0.34323	−0.3445	79	0.37531	0.37119	0.3753
52	-0.42972	-0.41933	−0.4314	48	0.4294	0.42768	0.4318
36	-0.40356	-0.37827	−0.4081	40	0.418656	0.41906	0.4224
32	-0.37977	-0.35344	−0.3839	36	0.407058	0.40917	0.4109
28	-0.34989	-0.32407	−0.3530	32	0.391374	0.3956	0.3949
1	0	0	0	1	0	0	0

Matlab Code URL

Link : https://github.com/RAAKASH/Intro-to-CFD-/tree/master/LDC
Check LDC - Post 3 for description of the Matlab code .

References

[1] Tamer A. AbdelMigid, Khalid M. Saqr, Mohamed A. Kotb, and Ahmed A. Aboelfarag. Revisiting the lid-driven cavity flow problem: Review and new steady state benchmarking results using {GPU} accelerated code. Alexandria Engineering Journal, 56(1):123 – 135, 2017.

[2] C. K. Aidun, N. G. Triantafillopoulos, and J. D. Benson. Global stability of a liddriven cavity with throughflow: Flow visualization studies. Physics of Fluids A: Fluid Dynamics, 3(9):2081–2091, 1991.

[3] O Botella and R Peyret. Benchmark spectral results on the lid-driven cavity flow. Computers & Fluids, 27(4):421–433, 1998.

[4] UKNG Ghia, Kirti N Ghia, and CT Shin. High-re solutions for incompressible flow using the navier-stokes equations and a multigrid method. Journal of computational physics, 48(3):387–411, 1982.

[5] K.A. Hoffman. Computational Fluid Dynamics for Engineers. Number v. 2. Engineering Education System, 1989.

Lid Driven Cavity - Mini Project (Post 1)

Lid Driven Cavity - Post 1

The following series of posts starting from this one gives you my approach in solving the Lid Driven Cavity , the numerical methods used and my collection of results up to Reynolds Number = 3200 . 

The series of posts also gives you my results in comparison with standard results that already exist for the problem such as the Ghia Ghia and Shin.  
This post gives you an outline of the approach used in solving the problem.
Before discussing the methods involved ,let me give you a brief introduction to what the problem actually is about. This problem is a standard benchmarking problem for testing any new software,code,etc. This consists of a rectangular or square cavity (2D or 3D) ,for now consider it to be a 2D square cavity whose top side is of infinite length driven with a constant velocity.

Image Courtesy : http://www.nacad.ufrj.br/~rnelias/gallery/cavity.html

Numerical Approach 

This problem has been solved using using the Stream Vorticity approach.

•  The Navier Stokes Equation has been solved using  Point Successive Over  Relaxation method.
• The Vorticity function has also been solved using the PSOR method.
• The Velocity functions have been solved by straight forward Central differencing scheme . 

  Algorithm used for Solving

1. Assuming a specific grid size  $N_x,N_y$  (best if they are equal).
2. Construct matrices for $\psi , \omega ,u,v$ using the above mentioned size using random values (best if values are zeros.)  
3. Assume appropriate boundary conditions for $\psi , \omega ,u,v$ . 
4. Solve for $\psi$ using the PSOR method using Vorticity function .
5.  Compute $u,v$ using central difference scheme .
6. Iterate the PSOR method for $\omega$  using Navier Stokes Equation (using the presolved $\psi,u,v$ .
   (Note : Don't solve for $\omega$ just run the PSOR over the entire grid just once.) 
7. Compute necessary Error Functions. Check for convergence.
8. If converged , STOP. Else go to Step 4.  

Note : In the above algorithm I haven't Explicitly mentioned what those Error functions are,there are many methods to check convergence . But I have mainly used four Parameters ,namely:

• Error 1 :RMSE $\left ( \frac{\psi  ^{k+1} _{i,j} -\psi  ^k _{i,j}}{\psi  ^k _{i,j}} \right )$ (where k is the iteration number)
  Variable name used in my code ("err").
• Error 2 :RMSE$\left ( \frac{\omega  ^{k+1} _{i,j} -\omega  ^k _{i,j}}{\omega  ^k _{i,j}}\right )$ .
  Variable name used in my code ("err1")
• Error 3 : Numerical solution accuracy of Vorticity equation .
  Variable name used in my code ("ERR1")
• Error 4 : Numerical solution accuracy of Navier stokes equation .
  Variable name used in my code ("ERR2")

It is necessary to note that out of these error parameters , Error 4 is the dominant parameter , this parameter would be used in monitoring the accuracy of the numerical scheme as well as convergence.

The following post would consist of the derivation of the stream vorticity functions , the math of the numerical schemes deployed. For now as a motivation for you to continue reading my blog I give you a beautiful graph of the vorticity contour solution for the $R_e = 600$.

Note: I still haven't introduced as of now the numerical schemes employed in Elliptic,hyperbolic PDEs (check timeline of blog),  which I soon hope to,but as of now I would explicitly mention the Numerical schemes used in the code.

1-D Parabolic Equation PDE schemes

1-D Heat conduction equation (with no internal heat generation)

$\frac{\delta T}{\delta t} = \alpha \left ( \frac{\delta ^2 T}{\delta x^2} \right )$

Before we dwell any further in the discretization of the above PDE we need to know 3 important terms.

1. Consistency
2. Stability 
3. Convergence

Here a finite difference scheme is said to be consistent when as $\Delta x,\Delta t \rightarrow 0$ , then the scheme  $\rightarrow$ to the original PDE.

A scheme is said to be stable when the ratio of the error propagation from one time step to another is less than $1$.

A scheme is said to converge when the above to conditions are met as well as the solution is close to the actual solution (analytic or numerical) of the PDE .

  

Note : Implicit schemes are mostly unconditionally stable ,while the explicit schemes are conditionally stable or unconditionally unstable. 

Explicit schemes

FTCS scheme

$\frac{T_{i} ^{n+1} -T_{i} ^{n}}{\Delta t} = \frac{T_{i+1} ^{n} -2T_{i} ^{n}+T_{i-1} ^{n}}{(\Delta x)^2}  + O(\Delta t,(\Delta x)^2)$

This scheme is conditionally stable , consistent ,convergent ,second order accurate in space,1st order accurate in time.
For more information on the stability criterion please check the tutorial 1 under the CFD tab.

CTCS scheme

$\frac{T_{i} ^{n+1} -T_{i} ^{n-1}}{2 \Delta t} = \frac{T_{i+1} ^{n} -2T_{i} ^{n}+T_{i-1} ^{n}}{(\Delta x)^2}  + O((\Delta t)^2,(\Delta x)^2)$

This scheme is unconditionally unstable. Hence this scheme is absolutely useless .Therefore the Dufort -Frankel scheme emerged.

Dufort Frankel scheme

$\frac{T_{i} ^{n+1} -T_{i} ^{n-1}}{2 \Delta t} = \frac{T_{i+1} ^{n} -2\left (\frac{T_{i} ^{n+1} +T_{i} ^{n-1}}{2}\right )+T_{i-1} ^{n}}{(\Delta x)^2}  + O((\Delta t)^2,(\Delta x)^2,(\frac{\Delta t}{\Delta x})^2)$

This scheme is unconditionally stable. Hence this scheme is brilliant . The error term can be easily derived , do this as an exercise.

Implicit schemes
BTCS scheme (laasonen scheme)

$\frac{T_{i} ^{n+1} -T_{i} ^{n}}{\Delta t} = \frac{T_{i+1} ^{n+1} -2T_{i} ^{n+1}+T_{i-1} ^{n+1}}{(\Delta x)^2}  + O(\Delta t,(\Delta x)^2)$

Unconditionally stable scheme . This scheme forms a tridiagonal matrix , which can be solved using matrix inversion ,or the Gauss elimination method or the TDMA (thomas) algorithm. For more information on the TDMA algorithm check tutorial 2 under the CFD tab .

Crank Nicolson scheme 

$\frac{T_{i} ^{n+1} -T_{i} ^{n}}{\Delta t} =  \frac{1}{2}\left (\frac{T_{i+1} ^{n+1} -2T_{i} ^{n+1}+T_{i-1} ^{n+1}}{(\Delta x)^2} +\frac{T_{i+1} ^{n} -2T_{i} ^{n}+T_{i-1} ^{n}}{(\Delta x)^2}    \right ) + O((\Delta t )^2,(\Delta x)^2)$

This is a really interesting scheme centered around the time step $n + 1/2$ . As can be seen the above scheme is second order accurate in time as well.
Unconditionally stable scheme . This scheme forms a tridiagonal matrix , which can be solved using matrix inversion ,or the Gauss elimination method or the TDMA (thomas) algorithm. 

Beta Formulation

Finally comes the $\beta $ formulation , this combines all the above mentioned schemes in a single shot .

$\frac{T_{i} ^{n+1} -T_{i} ^{n}}{\Delta t} =  \beta \left (\frac{T_{i+1} ^{n+1} -2T_{i} ^{n+1}+T_{i-1} ^{n+1}}{(\Delta x)^2}  \right )+(1-\beta) \left ( \frac{T_{i+1} ^{n} -2T_{i} ^{n}+T_{i-1} ^{n}}{(\Delta x)^2}   \right )  $

This scheme behaves like an FTCS scheme for $\beta < 1/2$
This scheme behaves like an BTCS scheme for $\beta > 1/2$
This scheme behaves like an Crank Nicolson scheme for $\beta = 1/2$ 

2D Parabolic PDEs

Here let us now consider a plate whose initial conditions are specified .
FTCS scheme
So Let us discretize the governing 2 -Dimensional Parabolic PDE using FTCS scheme .
$$\frac{T_{i,j}^{n+1}-T_{i,j}^{n}}{\Delta t} = \alpha \left (\frac{T_{i+1,j}^{n}-2T_{i,j}^{n}+T_{i-1,j}^{n}}{\Delta x ^2}  +\frac{T_{i,j+1}^{n}-2T_{i,j}^{n}+T_{i,j-1}^{n}}{\Delta y ^2} \right ) $$
Since this is an explicit scheme it becomes easy to calculate the solutions of $T$ as a function of time.(ie)

$$T_{i,j}^{n+1} = T_{i,j}^{n}+ \Delta t \alpha \left (\frac{T_{i+1,j}^{n}-2T_{i,j}^{n}+T_{i-1,j}^{n}}{\Delta x ^2}  +\frac{T_{i,j+1}^{n}-2T_{i,j}^{n}+T_{i,j-1}^{n}}{\Delta y ^2} \right )$$
For computational ease $\Delta x = \Delta y$ can be used. Now let us see the stability criterion of the scheme below.

Stability Condition

It can be derived to show that the stability for the 2 D heat conduction equation is similar to that of the 1-D heat equation ,hence 

$$ \alpha \Delta (t) \left (\frac{1}{(\Delta x )^2} + \frac{1}{(\Delta x )^2}  \right  )\leq \frac{1}{2}$$

Check Computational Fluid Dynamics by  Klaus A. Hoffmann (Author), Steve T. Chiang (Author) for reference.

This can be further reduced to 

$$ \alpha \Delta (t) \left  (\frac{1}{(\Delta x )^2} \right )\leq \frac{1}{4}$$ when $\Delta x = \Delta y$   

BTCS scheme

The same can be done in BTCS scheme to obtain 

$$\frac{T_{i,j}^{n+1}-T_{i,j}^{n}}{\Delta t} = \alpha \left (\frac{T_{i+1,j}^{n+1}-2T_{i,j}^{n+1}+T_{i-1,j}^{n+1}}{\Delta x ^2}  +\frac{T_{i,j+1}^{n+1}-2T_{i,j}^{n+1}+T_{i,j-1}^{n+1}}{\Delta y ^2} \right ) $$

The above equation as can be seen that it forms a pentadiagonal matrix which can be solved by matrix inversion which is computationally expensive ,hence this scheme wont be used most of the time. Hence emerges the ADI scheme

ADI scheme
The ADI (alternate direction implicit method) is the execution of the implicit schemes at one direction at time steps  $n,n +1/2$
$$\frac{T_{i,j}^{n+1/2}-T_{i,j}^{n}}{\Delta t} = \alpha \left (\frac{T_{i+1,j}^{n+1/2}-2T_{i,j}^{n+1/2}+T_{i-1,j}^{n+1/2}}{\Delta x ^2}  +\frac{T_{i,j+1}^{n}-2T_{i,j}^{n}+T_{i,j-1}^{n}}{\Delta y ^2} \right ) $$
$$\frac{T_{i,j}^{n+1}-T_{i,j}^{n+1/2}}{\Delta t} = \alpha \left (\frac{T_{i+1,j}^{n+1}-2T_{i,j}^{n+1}+T_{i-1,j}^{n+1}}{\Delta x ^2}  +\frac{T_{i,j+1}^{n+1/2}-2T_{i,j}^{n+1/2}+T_{i,j-1}^{n+1/2}}{\Delta y ^2} \right ) $$

The above equations as can be seen that it forms a tridiagonal matrix which can easily be solved by TDMA algorithm . The first equation is solved for finding the solutions for intermediate time step $n +1/2$ then this is used to find the solution at time step $n+1$ .