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...)

  • $\begingroup$ It seems that you are asking too many questions together. On one hand you are asking about the stability of your finite differences scheme, and on the other about fitting this to your experiment. $\endgroup$ – nicoguaro Sep 6 '17 at 20:24
  • $\begingroup$ Hi, the most important question is: "Is this a legitimate/viable approach to finding the diffusion constant?", basically if I make $\Delta t$ very small will it will converge to the physical situation so I can determine D from there by fitting the output to the experiment? $\endgroup$ – continuing_zeroes Sep 6 '17 at 21:24
  • 1
    $\begingroup$ You need to make both discretizations small enough for your solution to be close enough to the exact solution. Regarding your comment, if that differential equation models your experiment closely, then, in principle, the approach that you are proposing should work. $\endgroup$ – nicoguaro Sep 6 '17 at 22:24
  • 1
    $\begingroup$ "Is this issue a matter of convergence, and if so, how do I know how small to make Δx,Δt?" Yes. If the answer changes in an amount that is noticeable enough to you then it hasn't converged far enough. Of course, convergence is a limit to infinity, so you can never have it exact. But a good thing to do is to halve the timestep and see how different the solution is: this difference is a good estimate of the true error in some sense. You should keep halving until this difference is below your tolerance. Note that you should do this with both dx and dt. $\endgroup$ – Chris Rackauckas Sep 6 '17 at 22:56
  • 1
    $\begingroup$ Likely the easiest way to do this is via a method of lines. In this case you would discretize space by dx to get an ODE, and then give it to an ODE solver in say Julia or MATLAB. Adaptive ODE solvers have a tolerance, so you set that low enough so that you are comfortable with the time error. Then you just do the halving on dx until you're comfortable with the result. Adaptive solvers will naturally change the timestep so this is always stable, so this way you don't need to worry about the size of D affecting CFL. $\endgroup$ – Chris Rackauckas Sep 6 '17 at 22:58

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)

  • $\begingroup$ Thanks for the detailed explanations, definitely grateful to have more context to the problem! Confirmed that the code is consistent with the analytical solution and conditions described here, $e^{\frac{-x^2}{4Dt}} \frac{M}{4\pi D t}$. $\endgroup$ – continuing_zeroes Sep 7 '17 at 20:00

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. $$


Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.