Computational Physics: Finding the Diffusion Coefficient from the Discretized Diffusion Equation

Question

I'm pretty new to translating simulation to reality so please forgive the perhaps naive approach I'm taking here.

If we have a (quasi-2D) experimental video of a certain concentration changing with time, and we wish to determine the diffusion coefficient by comparing the video with a simulation, one can use the discretized diffusion equation (for a 2d square grid $\Delta x$ spacing) to map the behavior: $$\frac{\partial}{\partial t} \rho = D \nabla^2 \rho \Rightarrow \rho(i,j,t+\Delta t) = \rho(i,j,t) + \frac{D \Delta t}{\Delta x ^2}[\rho(i+\Delta x, j, t) + \rho(i-\Delta x, j, t) + \rho(i,j+\Delta x, t) + \rho(i, j-\Delta x, t) - 4\rho(i, j, t)]$$

So the parameter of interest is $\frac{D \Delta t}{\Delta x^2}$. I remember learning about the Courant number in Computational Physics, ($\frac{v \Delta t}{\Delta x}$) and it made sense the value couldn't be greater than 1, (or the real world was too fast for our simulation), but could be less than 1.

I've been given that the magical value for $\frac{D \Delta t}{\Delta x^2}$ is $\frac{1}{2}$ in order to maintain numerical stability.

Now what confuses me is that the only value that the simulation cares about is $\frac{D \Delta t}{\Delta x^2}$, but the value has 3 degrees of freedom (sorry for butchering that term).

So then, I'm allowed to choose $\Delta x$ and $\Delta t$, as long as $\frac{D \Delta t}{\Delta x^2} < .5$, correct?

But what confuses me is that if I take a given value of D and fix $\Delta x$, but change the value of $\Delta t$ by dividing $\frac{D \Delta t}{\Delta x^2}$ by 2 and increasing the number of loops through the diffusion process by 2, I do not receive the same outputs after n "seconds" have gone by (corresponding to x amount of relaxation steps for the simulation with $\frac{D \Delta t}{\Delta x^2}$ and 2x relaxation steps for the simulation with $\frac{1}{2} \frac{D \Delta t}{\Delta x^2}$. This seems to imply that if I were to attempt to find the real $D$ by fitting the experimental video to the simulation output, the value I determine for $D$ is dependent on $\Delta t$, which is certainly not good! And surely, $\Delta x$ would play a role as well then?

If anyone is kind enough to help me with any of these questions that would be really lovely!

1.) Is it down to an error in the code? (Or maybe the interpretation of it...)

2.) Is this issue a matter of convergence, and if so, how do I know how small to make $\Delta x, \Delta t$? Is there some proportionality with the test value for D?

3.) What's up with the $\frac{1}{2}$ value? (I've tried running values of $\frac{D \Delta t}{\Delta x^2}$ closer to $\frac{1}{2}$ and received errors...)

Answer 1

Your question contains many components, I will try to answer each individually

1.) Is it down to an error in the code? (Or maybe the interpretation of it...)

There are two main concepts in explicit finite difference solution of a linear parabolic equation (such as the transient diffusion equation). The first is stability. In order for your scheme to be stable, you need to have $$\frac{\Delta t D}{(\Delta x)^2} < 0.5 $$

This will you to a scheme that is stable in time and will not generate non-physical oscillations. The second one is accuracy. You are using first order in time explicit Euler scheme and second order in space centered discretization scheme. Therefore, the leading term of the norm (such as the $L^2$ norm) of your error will be : $$ e(\rho) \approx A_1O(\Delta x^2) + A_2 O (\Delta t)$$

Consequently, even if your scheme is stable, you are committing an error compared to the exact solution that is proportional to $\Delta x^2$ or $\Delta t$. Therefore, it is normal that your solution changes as you decrease $\Delta x$ and $\Delta t$ even if your stability condition is reached.

2) Gven a $\Delta x$, you can choose a $\Delta t$ for which your solution will not change significantly even if you decrease the time step by a factor of two. This means that your spatial discretization is now dominating the error term. If you do not find that you are converging to the same by solution by choosing smaller and smaller $\Delta t$ it means there is an error in your code. You can try to validate your code first using an analytical solution or the method of manufactured solution ( see this answer : In practice, what are the most useful ways to visualize 2d fluid flow, to tell what is happening in the simulation? Esp for verification and debugging )

3) The $\frac{1}{2}$ value arises from stability analysis. A crude explanation (not rigorous) can be given taking your equation in 1D:

$$\rho = D \nabla^2 \rho \Rightarrow \rho(i,t+\Delta t) = \rho(i,t) + \frac{D \Delta t}{\Delta x ^2}[\rho(i+\Delta x, t) + \rho(i-\Delta x, t) - 2\rho(i, t)]$$

If we take a simple example where $\rho(i,t)=0 \forall i \;\;except\;\; i=i_0$. Then the scheme for $i_0$ is :

$$ \rho(i_0,j,t+\Delta t) = \rho(i_0,t) + \frac{D \Delta t}{\Delta x ^2}[\rho(i_0+\Delta x, t) + \rho(i_0-\Delta x, t) - 2\rho(i_0, t)]$$

But given the initial conditions on $\rho$ we get : $$ \rho(i_0,j,t+\Delta t) = \rho(i_0,t) + \frac{D \Delta t}{\Delta x ^2}[ - 2\rho(i_0, t)]$$ This means that if $\frac{D \Delta t}{\Delta x ^2}>\frac{1}{2}$ then the concentration will become negative after a single iteration, something that is absolutely non-physical. Thus, for your scheme to remain stable you need $\frac{D \Delta t}{\Delta x ^2}<\frac{1}{2}$

This is really an hand waving explanation, but many books on the finite difference demonstrate much more rigourous ways to obtain this stability criterion (using the spectral radius of the finite difference scheme)

Answer 2

Just wanted to point out that if your grid spacing is $\Delta x$ and $\Delta y$ in $x$ and $y$ directions respectively, the stability criterion for the diffusion equation in 2-D is

$$ D \Delta t \left[\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} \right] \le \frac12. $$

So, when $\Delta x = \Delta y$, the above reduces to $$ \frac{D \Delta t}{\Delta x^2} \le \frac14. $$