I am reading this paper. In Example 2 from this paper, linear system of equation $Ax = b$ is given, where coefficient matrix $A$ has been generated by he five-point discretization of
the following second order elliptic partial differential equation:
$-\frac{\partial}{\partial x}(a\frac{\partial u}{\partial x})-\frac{\partial}{\partial y}(b\frac{\partial u}{\partial x}) + \frac{\partial}{\partial x}(cu)+ \frac{\partial}{\partial y}(du) +fu = 0 $$~~~~ (1)$
with $a(x, y)>0$, $b(x, y)>0$, $c(x, y)$, $d(x, y)$ and $f(x, y)$ defined on a unit square region $\omega = (0, 1)\times (0,1)$,
and Dirichlet boundary condition $u(x, y) = 0$ on boundary of $\omega$.
Now for the equation $(1)$ with $a(x, y) = b(x, y) = 1$, $c(x, y) = \cos(x/6)$, $d(x, y) = \sin(y/6)$, $f(x, y) = 1$ with three uniform meshes of $hx = hy = 1/11, 1/21, 1/31$ we get three matrices of orders $100\times 100$, $400\times 400$, $900\times 900$, where $hx$ and $hy$ refer to the mesh sizes in the $x$ direction
and $y$-direction, respectively.
I need help to generated matlab code for this problem. I have no idea how to apply matlab to generate these matrices. I would be very much grateful for any kind of help.
Thank you very much for your time.
Best Answer
You will have to augment whatever matrix you get as a result of applying a finite difference stencil discretizing the PDE with some equations representing the boundary conditions. The reason for this can be seen from looking at the one dimensional case. Say you had a just the ODE $u''(x) = 0$ on the domain $x \in (0,1)$. Say you want to discretize the domain into 5 intervals.
0 --- 0.2 --- 0.4 --- 0.6 --- 0.8 --- 1.0
Above I made a schematic diagram of the discretized domain. At the internal points you can get a finite difference approximation of the derivative as $\frac{u_{n-1}-2 u_n + u_{n+1}}{\Delta x}$ where $n$ is the index of the point the stencil is centered at. So you can get equations for the $x=0.2$ point where say $n=1$ and the $x= 0.4$ point where say $n=2$ up to the $x = 0.8$ point where $n=4$. But for $n=0$ and $n=5$ it won't work because the stencil would want to look $u_{-1}$ or $u_{6}$ which are outside our domain. So here is where the boundary conditions come in.
But for the internal points you can get 4 equations. In matrix form they would be
$$\left( \begin{array}{cccccc} 1 & -2 & 1 & 0 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 & 0 \\ 0 & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & 1 & -2 & 1 \\ \end{array} \right).\left( \begin{array}{c} u_0 \\ u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5 \\ \end{array} \right)=\left( \begin{array}{c} 0 \\ 0 \\ 0 \\ 0 \\ \end{array} \right)\Delta x$$
Adding the boundary conditions would make this a well posed problem. For example if $u(0) = 1$ and $u(1)=5$
$$\left( \begin{array}{cccccc} 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 & 0 \\ 0 & 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ \end{array} \right).\left( \begin{array}{c} u_0 \\ u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5 \\ \end{array} \right)=\left( \begin{array}{c} 1 \\ 0 \\ 0 \\ 0 \\ 0 \\ 5 \\ \end{array} \right)$$
It gets a little more complicated in two dimensions because you have to index the solutions values with a double index, but still list them in a one dimensional vector. Anyway if you posted an attempt you've made we could help you with it, but it's probaly unreasonable to ask people to write the matlab code for you from scratch.