# Implicit Finite difference scheme for a PDE with only one boundary

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.

• 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. Jan 27, 2014 at 16:10
• 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!
– DRG
Jan 27, 2014 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.

• Thanks for that; so essentially I keep altering my range and check that the solutions look similar over their length scales?
– DRG
Jan 27, 2014 at 17:10
• 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. Jan 27, 2014 at 17:49