# A Crank-Nicholson method and Robin boundary conditions

#### hunt_mat

Homework Helper
Summary
How does one implement the Robin condition for Crank-Nicholson?
I have the following PDE I wish to solve:
$$\frac{\partial u}{\partial t}=D\frac{\partial^{2}u}{\partial x^{2}}$$

With the following boundary conditions:

$$\frac{\partial u}{\partial x}(t,1)+u(t,1)=f(t),\quad u(t,0)=0$$

Now, I wish to do this via the Crank-Nicholson method and I would naively apply the method in the following way:
$$\frac{u_{i,N}-u_{i,N-1}}{\delta r}+u_{i,N}=f(t_{i})\equiv f_{i}$$

Now the stencil for the CN method is:
$$su_{i+1,j+1}-(2s+1)u_{i+1,j}+su_{i+1,j-1}=-(su_{i,j+1}-(2s+1)u_{i,j}+su_{i,j-1})$$

On the boundary this becomes:

$$su_{i+1,N}-(2s+1)u_{i+1,N-1}+su_{i+1,N-2}=-(su_{i,N}-(2s+1)u_{i,N-1}+su_{i,N-2})$$

I can write the boundary term as:

$$(1+\delta r)u_{i+1,N}=u_{i+1,N-1}+f_{i}\delta r$$

So I can then write:

$$\left(\frac{s}{1+\delta r}-(2s+1)\right)u_{i+1,N-1}+su_{i+1,N-2}=-(su_{i,N}-(2s+1)u_{i,N-1}+su_{i,N-2})-\frac{s\delta r}{1+\delta r}f_{i+1}$$

The system of equations is then:
$$A\mathbf{u}^{i+1}=-B\mathbf{u}^{i}=-\mathbf{b}$$

I've seen this done with a symmetric version of the first derivative, is that the right thing to do or is what I'm suggesting correct? I've got this coded up and it doesn't look the right sort of thing.

Related Differential Equations News on Phys.org

#### Chestermiller

Mentor
I would have used a central difference on the boundary condition, and combined it with the global difference equation at N to eliminate the value at N+1.

#### hunt_mat

Homework Helper
Yes, I realise why that is done, the method I initially proposed doesn't yield the derivative on the actual boundary. My code still doesn't seem to work as intended though.

#### mfig

You didn't supply the boundary condition at x=0, so I assume it is a constant. Incidentally, this looks like a Cartesian system given the PDE and x as a variable, but you have r in the boundary condition discretization so... I will assume the PDE is right.

If we use a central difference on the boundary then we end up with (I use superscripts for the time):

$\frac{u_{N+1}^{k} - u_{N-1}^{k}}{2\Delta x} + u_{N}^{k} = f_i$

This can be solved for the fake node at N+1 and put into the general formulation for the internal nodes. For 3 nodes, this gives the below system. Note for more nodes the second row is simply repeated with corresponding node adjustments in the u vector to form a tri-diagonal system.

$M = \begin{matrix} 1 & 0 & 0\\ -A & 2(1+A) & -A \\ 0&-2A & 2(1+AB) \end{matrix}$

$u^k = \begin{matrix} u_0^k \\ Au_0^k + 2(1-A)u_1^k + Au_2^k \\ 2Au_1^k + 2(1-AB)u_2^k+4A\Delta x f_2 \end{matrix}$

Where $A = \frac{D\Delta t}{(\Delta x)^2}$ and $B = 1 + \Delta x$. The updated values of u are given as solutions of $Mu^{k+1} = u^k$.

Last edited:

#### hunt_mat

Homework Helper
I thought I supplied the boundary condition at $x=0$, it was $u(t,0)=0$.

Last edited:

#### mfig

Yes, I see it now. I'm not sure how I missed it. At least it's a constant as I assumed.

#### hunt_mat

Homework Helper
I was hoping that my two methods of solution would a) work and b) give the same answer as one another.