I know the boundary condition is usually a tricky question. However, I am testing a finite-difference scheme for the equation of the form $$u_t=a(x)u_{xx}$$ that I know the analytical solution of. So if I specify the Dirichlet boundary condition, I can see the perfect convergence with a desired rate, but if I specify the value of the derivative I see the deterioration of accuracy of the numerical solution at the boundary(still have a good convergence in the middle of the grid). Both are exact from the analytical solution but looks like the derivative is not a good choice. What I think is that there should be only one solution that satisfies the initial condition and such a value of the "slope", or I am wrong? I might have a bug in my code but could not spot it so I wonder if imposing the partial derivative on the boundary instead of the value of the function could destroy uniqueness and this is why I am converging to a slightly different solution?
|
|
As Hui pointed out above, to apply Neumann boundary conditions correctly you should utilize ghost points and extend your $n^{th}$ order discretisation stencil to your domain boundaries. Utilising forward/backward difference or extrapolation at the boundary will degrade your solution. Assuming a $2^{nd}$ order central difference scheme the Neuman boundary condition, $u_x(t,x=0) = p_0(t)$, becomes (omitting $(t)$'s for convenience): $$ u_x\rvert_0 = p_0 \approx \frac{u_1 - u_{-1}}{2h} $$ $$ \Rightarrow u_{-1} = u_1 - 2hp_0 $$ Now, substitute for $u_{-1}$ into your discretised equation at the boundary, to obtain a Dirichelt BC $u_0$: $$ u_0(t) = a(0)\frac{u_1 - 2u_0 + u_{-1}}{h^2} $$ $$ \Rightarrow u_0 = a(0)\frac{2u_1 - 2hp_0}{h^2 + 2a(0)} $$ Applying the Neumann boundary condition in this way is consistent and is strictly equivalent to applying the analytic Dirichelt bc in the limit $h\rightarrow 0$. As for which is more accurate, for finite $h$, if you apply the Dirchlet condition at the boundary you impose exact agreement between the analytic and numerical solution at the boundary, but your estimation of the solutions derivative at this point will contain a $2^{nd}$ order error. Applying the Neumann boundary condition simply does the reverse, i.e. yields exact agreement between the analytic and numerical values of the derivative of your solution (estimated using your finite difference stencil) at the boundary, while estimations of $u$ at this point will contain a $2^{nd}$ order error. In the limit $h \rightarrow 0$ both methods should converge to the same unique solution. If you are not utilising ghost points then you are not applying the Neumann conditions in a consistent way and your result will indeed simply be less accurate then when you solve using the Dirichelt boundary conditions. |
|||||||
|
