Thread: Numerical solution of partial differential equation

1. Numerical solution of partial differential equation

Hi there. I am trying to self teach how to solve partial differential equations numerically using finite differences. I know this is a complex field, that requires much more knowledge of the theory than what I actually know, but anyway I wanted to try.

Anyway, I've tried to build my own differential equation, probably naively. I wanted to solve a first order partial differential equation, so I thought of taking the function:

$\displaystyle f(x,y)=sin(x-y)$

And then, taking derivatives I build the differential equation:

$\displaystyle \displaystyle \frac{\partial f(x,y)}{\partial x}+\frac{\partial f(x,y)}{\partial y}=0$, and I considered the interval $\displaystyle 0\leq x \leq 1, 0\leq y \leq 1$

Then I considered the boundary conditions $\displaystyle f(0,y)=\sin(-y)$, $\displaystyle f(1,y)=\sin(1-y)$,$\displaystyle f(x,0)=\sin(x)$ and $\displaystyle f(x,1)=\sin(x-1)$.

I think that everything is correct until here, but I'm not sure, if anything is wrong please tell me.

Then I considered centered differences (I needed for centered difference to have the previous and forward steps in the equation, in order to propagate the solution from the boundaries). I've considered equal spacings in the y and x directions, with $\displaystyle x_i=i\Delta x=ih$, $\displaystyle y_j=j\Delta y=j h$, and $\displaystyle h=\frac{1}{N+1}$, with N being the number of interior points in each axis, being $\displaystyle (x_0,y_j),(x_{N+1},y_j)$ and $\displaystyle (x_i,y_0),(x_i,y_{N+1})$ the points situated at the boundaries of the domain.

Under this conditions I've arrived to the set of equations:

$\displaystyle u_{i+1,j}-u_{i-1,j}+u_{i,j+1}-u_{i,j-1}=0$

Where I thought that $\displaystyle u_{i,j}\sim f(x_i,y_j)$

Then the discretized version of the boundary conditions would look like: [tex]u_{0,j}=\sin(-jh)[tex], $\displaystyle u_{N+1,j}=\sin(1-jh)$, $\displaystyle u_{i,0}=\sin(ih)$ and $\displaystyle u_{i,N+1}=\sin(ih-1)$.

But when I build the matrix and solve the linear system I get the result that the matrix is singular. So what would be the correct approach to this problem?

2. Re: Numerical solution of partial differential equation

Hey Ulysses.

You should look to see when the matrix you use to step through has a zero determinant.

If you are stepping through a partial differential equation, you will be doing something along the lines of a Taylor series expansion except you will be using linear operators [which are matrices] instead of normal derivatives.

Can you show us the scheme you are using?

As a hint - the Euler routine is y(x) = y(a) + y'(a)[x-a] + y''(a)[x-a]^2/2 + .... and instead of the derivatives you will be using linear operators and vectors instead of normal numbers.

3. Re: Numerical solution of partial differential equation

Hi, thank you very much for your kind response.

I am using, as I said, the same regular discretization in both directions with a step distance $\displaystyle h=\frac{1}{N+1}$: $\displaystyle x_i=i\Delta x=ih$ and $\displaystyle y_j=j\Delta y=jh$

with $\displaystyle 0\leq x,y \leq 1$, so that: $\displaystyle (x_0,y_0)=(0,0)$ and $\displaystyle (x_{N+1},y_{N+1})=(1,1)$

In the discrete case I have this boundary conditions:
$\displaystyle u_{0,j}=\sin(-y_j)=\sin(-jh)$, $\displaystyle u_{N+1,j}=\sin(1-y_j)=\sin(1-jh)$ for all $\displaystyle j=0,1,2,...,N+1$
$\displaystyle u_{i,0}=\sin(x_i)=\sin(ih)$, $\displaystyle u_{i,N+1}=\sin(x_i-1)=\sin(1-ih)$ for all $\displaystyle i=0,1,2,...,N+1$

From centered finite differences, the differential equation is:

$\displaystyle u_{i+1,j}-u_{i-1,j}+u_{i,j+1}-u_{i,j-1}=0$ where I considered this equation for all i,j being interior points on the domain.

for the case N=2 from this set of equations I get this linear system:

$\displaystyle \begin{bmatrix} 0 & 1 & 1 & 0 \\ -1 & 0 & 0 & 1 \\ -1 & 0 & 0 & 1 \\ 0 & -1 & -1 & 0 \end{bmatrix} \left[ \begin{array}{c} u_{1,1} \\ u_{2,1} \\ u_{1,2} \\ u_{2,2} \end{array} \right] = \left[ \begin{array}{c} u_{0,1}+u_{1,0} \\ u_{2,0}+u_{3,1} \\ u_{0,2}+u_{1,3} \\ -u_{2,3}-u_{3,2} \end{array} \right]$

I have actually tried to solve the linear system for the N=3 case, which contains a 9x9 matrix, but I think that if the problem is in the algorithm, the case N=2 should equally serve to discuss the problem.