I have the following parabolic equation:
$$ \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} + \frac{\partial ^2u}{\partial y^2} $$
over domain $(x,y)\in [0,10] \times [0,10]$ where $\Delta x = \Delta y = 1 \times 10^{-2}$ and $\Delta t = 2.4 \times 10^{-3}$. Initial and boundary conditions are:
$$ \left\{ \begin{array}{ll} u(x,y,0) = 10\\ u(0,y,t) = 4\\ u(10,y,t) = 4\\ u(x,0,t) = 9\\ u(x,10,t) = 9\\ \end{array} \right. $$ I tried to rewrite the equation as follows:
$$ \frac{u_{i,j}^{k+1}-u_{i,j}^{k}}{\Delta t} = \frac{u_{i+1,j}^{k+1}-2u_{i,j}^{k+1}+u_{i-1,j}^{k}}{\Delta x^2} + \frac{u_{i,j+1}^{k+1}-2u_{i,j}^{k+1}+u_{i,j-1}^{k}}{\Delta y^2} $$
I initially wanted to write the equation as a tridiagonal matrix:
$$ -\frac{1}{\Delta x^2}(u_{i+1,j}^{k+1}+u_{i-1,j}^{k+1}) - \frac{1}{\Delta y^2} (u_{i,j+1}^{k+1}+u_{i,j-1}^{k+1})+(\frac{1}{\Delta t}+\frac{2}{\Delta x^2}+2\frac{1}{\Delta y^2})u_{i,j}^{k+1}= \frac{u_{i,j}^{k}}{\Delta t}. $$
However, I didn't found a way to do so. I am not sure how to solve this, but is it possible that I have to solve this by iteration like:
$$ u_{i,j}^{k+1}= \Big(\frac{u_{i,j}^{k}}{\Delta t}+\frac{1}{\Delta x^2}(u_{i+1,j}^{k+1}+u_{i-1,j}^{k+1}) + \frac{1}{\Delta y^2} (u_{i,j+1}^{k+1}+u_{i,j-1}^{k+1})\Big)(\frac{1}{\Delta t}+\frac{2}{\Delta x^2}+2\frac{1}{\Delta y^2})^{-1}. $$
As user_of_math suggested, you should be consistent with your $k$ index.Let $k$ denote the index for time, $i$ for $x$-direction in space and $j$ for $y$-direction in space. Then you use the wrong discretization. It should always be $k+1$ as the time index except for the time derivative since your scheme should be implicit. \begin{align} \frac{u^{k+1}_{i,j}-u^{k}_{i,j}}{\Delta t} = \frac{u^{k+1}_{i+1,j}-2u^{k+1}_{i,j}+u^{k+1}_{i-1,j}}{\Delta x^2} + \frac{u^{k+1}_{i,j+1}-2u^{k+1}_{i,j}+u^{k+1}_{i,j-1}}{\Delta y^2} \end{align} or when you want to solve for the next time step \begin{align} u^{k+1}_{i,j} = u^{k}_{i,j}+\Delta t\Bigl(\frac{u^{k+1}_{i+1,j}-2u^{k+1}_{i,j}+u^{k}_{i-1,j}}{\Delta x^2} + \frac{u^{k+1}_{i,j+1}-2u^{k+1}_{i,j}+u^{k+1}_{i,j-1}}{\Delta x^2}\Bigr) \end{align} and since $\Delta x =\Delta y$ \begin{align} &u^{k+1}_{i,j} = u^{k}_{i,j}+\frac{\Delta t}{\Delta x^2}\Bigl({u^{k+1}_{i+1,j}+u^{k+1}_{i-1,j}} -4u^{k+1}_{i,j} + {u^{k+1}_{i,j+1}+u^{k+1}_{i,j-1}}\Bigr) \\ \Leftrightarrow& u^{k+1}_{i,j}-\frac{\Delta t}{\Delta x^2}\Bigl({u^{k+1}_{i+1,j}+u^{k+1}_{i-1,j}} -4u^{k+1}_{i,j} + {u^{k+1}_{i,j+1}+u^{k+1}_{i,j-1}}\Bigr) =u^k_{i,j} \end{align} You get this for each pair $(i,j)$ giving you $N_x \times N_y$ equations. Assume $N=N_x=N_y$. Now you need to solve a linear system in each time step. You have to come up with some ordering for your matrix and right hand side. The standard ordering would give you a matrix $A \in \mathbb{R}^{N^2\times N^2}$ which looks like this: \begin{align} A = \begin{pmatrix} \tilde A &-\tilde I &0 & ...&...&0 \\ -\tilde I & \tilde A&-\tilde I &0 & .. & 0 \\ 0 &-\tilde I & \tilde A &-\tilde I & 0 &..\\ &&& \ddots \\ 0 &...&...&0&-\tilde I & \tilde A \end{pmatrix} \end{align} Where $\tilde I$ is the $n\times n$ identiy matrix, scaled by $c:=\Delta t / \Delta x^2$ and $\tilde A$ is an $n\times n$ matrix. \begin{align} \tilde A = \begin{pmatrix} 1+c&-c &0 & ...&...&0 \\ -c & 1+c&-c &0 & .. & 0 \\ 0 &-c & 1+c &-c & 0 &..\\ &&& \ddots \\ 0 &...&...&0&-c & 1+c \end{pmatrix} \end{align} Now you solve \begin{align} Au^{k+1}=u^k \end{align} during each time step.
Edit: Of course you have to take care of the boundary values during the computation. I did not consider this.