Forum Navigation: Select WHO'S ON: 02:22 PM General Forum Technical Forum Economics Forum Numerical Methods Forum Trading Forum The Quantitative Finance FAQs Proje... Student Forum Book And Research Paper Forum Programming and Software Forum The Quantitative Finance Code Libra... Careers Forum Events Board Brainteaser Forum Off Topic Forum and Website Bugs and Suggesti... Wilmott / BookMap Competition

 new topic
 search
 jobs board
 magazine
 news
 help
 file share
 audio visual
 articles
 blogs
 wiki
 forums
 home

 FORUMS > Numerical Methods Forum < refresh >
 Topic Title: Numerical solution of the Fokker-Plank equation Created On Wed Aug 11, 10 08:23 PM Topic View: Branch View Threaded (All Messages) Threaded (Single Messages) Linear

 Pages: [ 1 2 >> Next ]

MJPL
Junior Member

Posts: 7
Joined: Aug 2010

Wed Aug 11, 10 08:23 PM

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?

M

antiser
Member

Posts: 34
Joined: Jun 2009

Wed Aug 11, 10 08:52 PM

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.

vandervolt
Member

Posts: 51
Joined: Apr 2007

Thu Aug 12, 10 06:45 AM

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

MJPL
Junior Member

Posts: 7
Joined: Aug 2010

Thu Aug 12, 10 09:13 AM

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

SgtPepper
Member

Posts: 24
Joined: Nov 2004

Thu Aug 12, 10 01:19 PM

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

MJPL
Junior Member

Posts: 7
Joined: Aug 2010

Thu Aug 12, 10 01:25 PM

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

SgtPepper
Member

Posts: 24
Joined: Nov 2004

Thu Aug 12, 10 02:01 PM

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.

Orbit
Senior Member

Posts: 471
Joined: Oct 2003

Thu Aug 12, 10 02:54 PM

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.

Alan
Senior Member

Posts: 8428
Joined: Dec 2001

Thu Aug 12, 10 03:09 PM

Quote

Originally posted by: MJPL
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...

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 problems
can be approached from either.

My suggestions:

-First, categorize the boundary behavior using
the Feller criteria as explained in Karlin & Taylor or other standard sources.

-Second, develop the exact soln behavior as the boundaries are approached
for both the backwards and forward (FP) eqn. This may suggest other
boundaries 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 mesh
pts to the right of p as to the left. Try it on both the forward and backward problems., perhaps
with 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, and
your favorite differencing scheme. For example, try y = x/(x+p) which
has very nice properties.

-------------------------
www.financepress.com

Edited: Thu Aug 12, 10 at 04:44 PM by Alan

vandervolt
Member

Posts: 51
Joined: Apr 2007

Fri Aug 13, 10 04:08 AM

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 is

dx = 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_x

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

MJPL
Junior Member

Posts: 7
Joined: Aug 2010

Fri Aug 13, 10 01:32 PM

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

Cuchulainn
Senior Member

Posts: 47582
Joined: Jul 2004

Fri Aug 13, 10 03:58 PM

Quote

Originally posted by: MJPL
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?

M

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

-------------------------
www.datasimfinancial.com

Edited: Fri Aug 13, 10 at 04:05 PM by Cuchulainn

MJPL
Junior Member

Posts: 7
Joined: Aug 2010

Fri Aug 13, 10 04:35 PM

Quote

Originally posted by: Cuchulainn

Fitting 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)?

Quote

I 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?

Cuchulainn
Senior Member

Posts: 47582
Joined: Jul 2004

Sat Aug 14, 10 04:55 PM

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?

-------------------------
www.datasimfinancial.com

MJPL
Junior Member

Posts: 7
Joined: Aug 2010

Sat Aug 14, 10 05:28 PM

Quote

Originally posted by: Cuchulainn
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?

No, I don't actually. Can you suggest a few?

Thanks so much.... M

Cuchulainn
Senior Member

Posts: 47582
Joined: Jul 2004

Tue Aug 17, 10 01:41 PM

Quote

Originally posted by: MJPL
Quote

Originally posted by: Cuchulainn
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?

No, I don't actually. Can you suggest a few?

Thanks so much.... M

Thinking out loud. Since FPE is 2nd order you need 2 extra conditions, so where will they come from? E.g. at x = 0 we could say the pde is satisfied there or Dirichlet/Neumann is defined there. That leaves 1 condition which could be the intergal == 0 as you stated.

It could be some variational calculus problem See Gelfand and Fomin.

You could try a fd scheme and show it has a solution and then let h ->0 to prove the pde is well-posed.

I am sure this pde has been posed somewhere.

-------------------------
www.datasimfinancial.com

Alan
Senior Member

Posts: 8428
Joined: Dec 2001

Tue Aug 17, 10 02:02 PM

Quote

Originally posted by: Cuchulainn
I am sure this pde has been posed somewhere.

It is a running example in Karlin & Taylor's 'Second Course ...' -- the Wright-Fisher gene frequency diffusion model.
Googling Wright-Fisher diffusion turns up about 16 pages of refs.

-------------------------
www.financepress.com

Edited: Tue Aug 17, 10 at 02:07 PM by Alan

MJPL
Junior Member

Posts: 7
Joined: Aug 2010

Tue Aug 17, 10 02:16 PM

Quote

Originally posted by: Alan
It is a running example in Karlin & Taylor's 'Second Course ...' -- the Wright-Fisher gene frequency diffusion model.
Googling Wright-Fisher diffusion turns up about 16 pages of refs.

I am very familiar with Karlin and Taylor's book and with the mathematical theory of the problem. I am just struggling a bit with the numerical solution (not covered in Karlin and Taylor!).

Edited: Tue Aug 17, 10 at 02:18 PM by MJPL

Alan
Senior Member

Posts: 8428
Joined: Dec 2001

Tue Aug 17, 10 03:55 PM

try y = x/(x+p)

-------------------------
www.financepress.com

Cuchulainn
Senior Member

Posts: 47582
Joined: Jul 2004

Wed Aug 18, 10 01:34 PM

Quote

Originally posted by: Alan
Quote

Originally posted by: Cuchulainn
I am sure this pde has been posed somewhere.

It is a running example in Karlin & Taylor's 'Second Course ...' -- the Wright-Fisher gene frequency diffusion model.
Googling Wright-Fisher diffusion turns up about 16 pages of refs.

nice. Thanks. Knowing the problem helps finding a solution

-------------------------
www.datasimfinancial.com