Originally asked this on math.stackexchange, but I figure it's also appropriate here. I'm reading through some finite difference code for a diffusion equation and came across something odd for the boundary conditions, and I was wondering about the validity of the method. I'm familiar with the finite difference method, but I've never seen this certain method before.

The equation I am solving is

$$ \frac{\partial f}{\partial t}=C(x)\frac{\partial^2 f}{\partial x^2}$$

with boundary condition

$$ \frac{\partial f}{\partial x}=0\ \text{on}\ x=0,1 $$

across time, $t$, and position, $x$, with position-dependent diffusion constant, $C(x)$. It's being solved in an explicit way using the forward difference for time and the second-order central difference for space

$$ \frac{f^\text{new}_i-f^\text{old}_i}{\Delta t}=C_i\frac{f^\text{old}_{i-1}-2f^\text{old}_{i}+f^\text{old}_{i+1}}{(\Delta x)^2}$$

along discretised space $\{x_1,\dots,x_n\}$. Rearranging the above equation gives

\begin{align} f^\text{new}_i=&f^\text{old}_i + \frac{C_i\Delta t}{(\Delta x)^2}\Big(f^\text{old}_{i-1}-2f^\text{old}_{i}+f^\text{old}_{i+1}\Big)\\ =&f^\text{old}_{i-1}(S_i)+f^\text{old}_{i}(1-2S_i)+f^\text{old}_{i+1}(S_i) \end{align}

where $S_i=\frac{C_i\Delta t}{(\Delta x)^2}$.

This is simply represented in the code. For the boundary at $x_1$ (and at $x_n$, but I'll only talk about the first as they are similar), there are some strange things going on. Usually I would do a Taylor series expansion of function $f$ at points $x_2$ and $x_3$

$$ f_2=f_1+\frac{\partial f}{\partial x}\bigg\rvert_i \Delta x + \frac{1}{2} \frac{\partial^2 f}{\partial x^2}\bigg\rvert_i (\Delta x)^2 + \mathcal{O}((\Delta x)^3)$$ $$ f_3=f_1+\frac{\partial f}{\partial x}\bigg\rvert_i (2\Delta x) + \frac{1}{2} \frac{\partial^2 f}{\partial x^2}\bigg\rvert_i (2\Delta x)^2 + \mathcal{O}((\Delta x)^3)$$

which, due to boundary condition $\frac{\partial f}{\partial x}=0$ can be rearranged to

$$ \frac{\partial^2 f}{\partial x^2}\bigg\rvert_i\approx\frac{2}{(\Delta x)^2}\Big(f_2 - f_1\Big)\tag{$A$}$$ $$ \frac{\partial^2 f}{\partial x^2}\bigg\rvert_i\approx\frac{1}{2(\Delta x)^2}\Big(f_3 - f_1\Big)\tag{$B$}.$$

Usually, I'd average $A$ and $B$ to determine how the boundary conditions are applied in the eventual matrix solution

$$ \frac{\partial^2 f}{\partial x^2}\bigg\rvert_i\approx\frac{1}{(\Delta x)^2}\Big(-\frac{5}{4}f_1 + f_2 + \frac{1}{4}f_3 \Big).$$

However, the code that I'm reading does something different and takes the boundary as a weighted sum of $A/5+4B/5$, as shown below

$$\frac{A+4B}{5} \Rightarrow \frac{\partial^2 f}{\partial x^2}\bigg\rvert_i\approx\frac{1}{(\Delta x)^2}\Big(-\frac{4}{5}f_1 + \frac{2}{5}f_2 + \frac{2}{5}f_3 \Big).$$

I don't think this as being inherently wrong, as the LHS is still $\frac{\partial^2 f}{\partial x^2}\big\rvert_i$, but I'm not sure. Is this a valid step to make? Is it ok to "weight" results from different nodes, and if so, in what situations would weighting or not weighting be useful?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.