Dear all,I need to obtain an accurate numerical solution to the Fokker-Plank/ forward equation:for , subject to and zero-flux conditions on the space boundaries. I am specifically interested in the case where p is very close to 0 or 1. I have attempted to solve this using both the Crank-Nicolson method and the exponentially-fitted finite difference method of Duffy (with uniform and non-uniform meshes). Unfortunately neither method produces accurate estimates when is concentrated near the boundaries. Does anyone know of superior numerical routines for solving this equation when the initial condition is concentrated near the boundaries?Thanks in advanceM

Finite element method would be great for that as you can mesh the point of interest with needed accuracy.MATLAB pde toolbox should be able to do this, although it's been a while since I used it.I simulate similar problems in different field using COMSOL.For 1D as you have you should be able to code it rather quickly.Look for references in wiki page.

It is a little difficult to say without knowing what your a(x) and b(x) terms are. That said, I currently use exponential fitting and Rannacher time stepping for a Fokker-Planck equation with a = a(x,t) and b= b(x) with two zero flux boundaries and it works nicely. you may get nasty behaviour about the boundaries depending on a(x) and b(x). perhaps in these cases you need very fine grids which are impractical in Matlab (systems of equation solver to slow...). also, as you mention Duffy's work be careful of the delta function IC with crank-nicolson.. things can go belly-up...also, antiser is right - do you have access to Matlab? the PDE solver pdepe is (imo) solid and easy to use. Just write your PDE in conservative form and away you go. What are your a(x) and b(x)?

Thanks Vandervolt and Antiser!I haven't tried the finite element method before and will try and get Matlab from someone in the department to give it a shot. @ Vandervolt: In my application (which is actually a population genetics problem so I left out the details in this forum!) I have a(x) = c x(1-x) and b(x) = x(1-x), where c is a finite (positive or negative) constant. Is there a general rule as to when the values of a(x) and b(x) would produce bad boundary behaviour?I have tried finer uniform grids and non-uniform grids which are finer near the boundaries but that doesn't seem to help. I am aware of the problems with CN when the IC is a delta function which is why I swapped to exponential fitting. Let me have a go with Matlab and see it that works...

Hi,could you please detail what you mean with "zero-flux conditions on the space boundaries" ? Does it mean that for x = 0 and for x = 1 ? Is there any justification ? Many thanks.SgtPepper

Hi SgtPepper,That is precisely what I mean. The function is a density and this condition ensures that it always integrates to one. Thanks,M

Thanks MJPL,please allow me to insist, but what should I take into account between a condition on the Spot such as and a condition on x = log(S) such as ? Thanks again.

That's the problem, such conditions produce unstable oscillations for small x but even moreso for x when tau is almost zero. This can be very frustrating but as the previous posters have indicated I think FEM can help.

QuoteOriginally posted by: MJPLThanks Vandervolt and Antiser!I haven't tried the finite element method before and will try and get Matlab from someone in the department to give it a shot. @ Vandervolt: In my application (which is actually a population genetics problem so I left out the details in this forum!) I have a(x) = c x(1-x) and b(x) = x(1-x), where c is a finite (positive or negative) constant. Is there a general rule as to when the values of a(x) and b(x) would produce bad boundary behaviour?I have tried finer uniform grids and non-uniform grids which are finer near the boundaries but that doesn't seem to help. I am aware of the problems with CN when the IC is a delta function which is why I swapped to exponential fitting. Let me have a go with Matlab and see it that works...Your boundaries are singular. If the particle starts there, it cannot escape. On the other hand, if it starts in the interior, the boundaries may be unreachable.A zero-flux condition may be illegal (i.e. flat out wrong) or simply not the best choice.You may or may not need to truncate your boundaries to (eps,1-eps), where eps is a small number.Also, the backwards Kolmogorov eqn may be a better choice, as most problemscan be approached from either.My suggestions:-First, categorize the boundary behavior usingthe Feller criteria as explained in Karlin & Taylor or other standard sources. -Second, develop the exact soln behavior as the boundaries are approachedfor both the backwards and forward (FP) eqn. This may suggest otherboundaries conditions, such as simple absorption, or simply"no boundary condition" (i.e. just apply the pde there like a reg. mesh pt)-Try this mesh with p small: 2-sided uniform with (approximately) the same number of meshpts to the right of p as to the left. Try it on both the forward and backward problems., perhapswith boundary truncation. A somewhat similar (and really cleaner) alternative to this is to map to a new variable y = y(x) = y(x;p). Here y ranges from 0 to c, where c is close to 1, but the "hotspot" x=p is mapped to y=0.5. Then, just use a uniform mesh in y, andyour favorite differencing scheme. For example, try y = x/(x+p) whichhas very nice properties.

I have a related question regarding the boundary conditions MJPL is using...If you want "zero-flux" boundaries, I would have set not \phi_x = 0, where \phi_x is the partial derivative of \phi wrt x. Although it does make sense that if you set \phi_x = 0 on the boundaries total probability is still being conserved. [and to check I take a simple case where b(x) = b and a(x) = x(1-x) and numerically probability is conserved in both cases...]So the question is - which is the "correct" way to implement a zero-flux boundary? More importantly, the underlying process in MJPL's case isdx = x(1-x)dt + \sqrt[ x(1-x) ] dW.So what behaviour does this SDE have in the two cases (i) Zero flux (ii) \phi_xProbability is being conserved in both cases, but since the boundary conditions are different, surely the boundary behaviour is different?Thanks MJPL for bringing up the question, very interesting and now you've got me confused... btw Im thinking of a setting where the left boundary is at x=0, and the right boundary at any x>0 (not necessarily x=1).

Thanks to everyone for their replies! Your comments have been very helpful. One of the suggestions was that I solve the pde using the finite element method implemented in the pde toolbox of Matlab. However, I am not sure how to specify this form in the "PDE specification" menu (this is the first time I am using Matlab!). Are there any Matlab gurus out there that know how?Thanks,M

QuoteOriginally posted by: MJPLDear all,I need to obtain an accurate numerical solution to the Fokker-Plank/ forward equation:for , subject to and zero-flux conditions on the space boundaries. I am specifically interested in the case where p is very close to 0 or 1. I have attempted to solve this using both the Crank-Nicolson method and the exponentially-fitted finite difference method of Duffy (with uniform and non-uniform meshes). Unfortunately neither method produces accurate estimates when is concentrated near the boundaries. Does anyone know of superior numerical routines for solving this equation when the initial condition is concentrated near the boundaries?Thanks in advanceMFitting won't help here because it is not a convection-dominated problem. Nor the CN..I am wondering why x is in the interval [0,1] Normally it is on the real line or the half real line. If your Pde is on the real line (no pde Boundary conditions), one approach is to transform x to y = tanh(x) in interval [-1,1] like I do in my ADE article and then BC become easier. I have done thus for other pdes but not FPE yet. I would be more inclined to give Dirichlet BC..If you take the heat pde you see solution --> 0 at infinities, not zero flux BC which is mathematically incorrect. It has an exact solution. I would not be a fan of FEM for this. It is a sledgeHammer.//!!The basic assumption is wrong; Pdes on real line do not have/need BC, only the numerical approximations need BC.

QuoteOriginally posted by: CuchulainnFitting won't help here because it is not a convection-dominated problem. Nor the CN..Surely this would depend on the relative magnitudes of a(x) and b(x)?QuoteI am wondering why x is in the interval [0,1] Normally it is on the real line or the half real line. Perhaps some context would help here. I am modelling the frequency of a mutation in a population - so x is a proportion bound between 0 and 1. The probability of the population having a frequency x at time t is described by the density function:where p_0 = Pr[X=0] and p_1 = Pr[X=1]. The function does not have any singularities and evolves according to the Fokker-Planck equation subject to boundary conditions that ensure that for all t.Any further suggestions on how to obtain an accurate numerical solution for in this context?

the only suggestion at this stage is that the problem looks like a variational problem with an integrability constraint. ??Do you have references to similar problems?

QuoteOriginally posted by: Cuchulainnthe only suggestion at this stage is that the problem looks like a variational problem with an integrability constraint. ??Do you have references to similar problems?No, I don't actually. Can you suggest a few?Thanks so much.... M