I am looking at a few reaction-diffusion equations of the form

$\frac{dP}{dt} = D\left(\frac{d^2P}{dr^2} + \frac{2}{r}\frac{dP}{dr}\right) - a(P)$

I know the initial conditions and the boundary value at one end. I also know the function steadily falls and eventually hits zero at the second unknown boundary. To model this behaviour and indirectly find this unknown boundary, I coded a solver using an explicit finite difference scheme, with an added condition in the loop that changed forced any negative values to zero. This gave me the correct result, but due to the CFL condition I've had to use a tiny time step ($\Delta t = 0.001$) when I'm more interested in the system's behaviour over several hours or even days.

I was looking into using implicit finite difference methods or CN methods so I may increase the time step, but my (limited) understanding of this implies I'd have to solve a system of equations which would include the unknown boundary, which I will not know exactly. Is it possible to work around this, or will implicit methods fail? If it is possible to work around it, could anyone suggest a good method and how I would implement this? Thanks in advance, I'm quite new to numerical methods and would appreciate any suggestions.

  • $\begingroup$ Could you explain more what you did with the explicit scheme with the second boundary? I don't see how what you currently said would in any way handle the boundary. $\endgroup$ – Godric Seer Jan 27 '14 at 16:10
  • $\begingroup$ Sure; First I made a grid. I know initially all points on the grid except the boundary values are 0, so I set them to that. I knew at P(x0,t) = po , so I specified that. I then used FD method over 500 spatial steps and 600 time steps (10 minutes) to compute P(i,j +1) ; I used an IF loop to check P(i,j+1) was positive or zero; if it wasn't, it was set to zero before the loop incremented. This method worked quite well, but I suspect I can't make an implicit analog! $\endgroup$ – DRG Jan 27 '14 at 16:25

I am not entirely certain that your method works correctly. The answers may look correct but unless you have an analytic solution to check against you could be introducing errors. Also, how do you handle the 2nd derivative at the far boundary?

The boundary condition you have is:

$ \lim \limits_{x \to \infty} P = 0$

with a domain of $[x_0, \infty)$. A more common way to handle this numerically is to truncate your domain to $[x_0, x_1]$ and hold the boundary condition that

$ P(x_1) = 0$

This works fine, except for you don't know if your truncated domain has affected your answer. So then you resolve it on $[x_0, 2 x_1]$ (assuming $x_1 >> x_0$ and this about doubles your domain). If your answer in the areas of interest are close, then you know you are done.

If you make this change, then switching to an implicit solver would be trivial, since you have eliminated the if condition inside your solver.

  • $\begingroup$ Thanks for that; so essentially I keep altering my range and check that the solutions look similar over their length scales? $\endgroup$ – DRG Jan 27 '14 at 17:10
  • $\begingroup$ Yes. You know it is 0 out somewhere, so you just apply it somewhere and check that it is far enough out to not affect your solution. $\endgroup$ – Godric Seer Jan 27 '14 at 17:49

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.