Crank-Nicholson method and Robin boundary conditions

[hunt_mat]

TL;DR

How does one implement the Robin condition for Crank-Nicholson?

I have the following PDE I wish to solve:
<br /> \frac{\partial u}{\partial t}=D\frac{\partial^{2}u}{\partial x^{2}}<br />

With the following boundary conditions:

<br /> \frac{\partial u}{\partial x}(t,1)+u(t,1)=f(t),\quad u(t,0)=0<br />

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

Now the stencil for the CN method is:
<br /> 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})<br />

On the boundary this becomes:

<br /> 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})<br />

I can write the boundary term as:

<br /> (1+\delta r)u_{i+1,N}=u_{i+1,N-1}+f_{i}\delta r<br />

So I can then write:

<br /> \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}<br />

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

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.

 

[Chestermiller]

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]

Yes, I realize 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$.

 

[hunt_mat]

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

 

[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]

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

 

[adomanim1]

with solution A1(r,ϵ)A1(r,ϵ). In particular, for ϵ=0ϵ=0 you get the solution A1(r,0)A1(r,0) which is a solution for ddrA=F(A,r)ddrA=F(A,r), so you should get the same solutions (if they agree at a point and the differential equation behaves nicely enough to get a unique solution).Tutuapp[/size][/color] 9apps[/size][/color] ShowBox[/size][/color]

 

[hunt_mat]

This has now been solved and is a dead thread.