All New Wilmott Jobs Board                     (g)

Forum Navigation:

magazine

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:

Pages: [ 1 2 >> Next ]
View thread in raw text format


MJPL
Junior Member

Posts: 7
Joined: Aug 2010

Wed Aug 11, 10 08:23 PM
User is offline

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 advance

M
 
Reply
   
Quote
   
Top
   
Bottom
     



antiser
Member

Posts: 34
Joined: Jun 2009

Wed Aug 11, 10 08:52 PM
User is offline

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.
 
Reply
   
Quote
   
Top
   
Bottom
     



vandervolt
Member

Posts: 51
Joined: Apr 2007

Thu Aug 12, 10 06:45 AM
User is offline

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)?
 
Reply
   
Quote
   
Top
   
Bottom
     



MJPL
Junior Member

Posts: 7
Joined: Aug 2010

Thu Aug 12, 10 09:13 AM
User is offline

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...
 
Reply
   
Quote
   
Top
   
Bottom
     



SgtPepper
Member

Posts: 24
Joined: Nov 2004

Thu Aug 12, 10 01:19 PM
User is offline

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
 
Reply
   
Quote
   
Top
   
Bottom
     



MJPL
Junior Member

Posts: 7
Joined: Aug 2010

Thu Aug 12, 10 01:25 PM
User is offline

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
 
Reply
   
Quote
   
Top
   
Bottom
     



SgtPepper
Member

Posts: 24
Joined: Nov 2004

Thu Aug 12, 10 02:01 PM
User is offline

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.
 
Reply
   
Quote
   
Top
   
Bottom
     



Orbit
Senior Member

Posts: 277
Joined: Oct 2003

Thu Aug 12, 10 02:54 PM
User is offline

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.
 
Reply
   
Quote
   
Top
   
Bottom
     



Alan
Senior Member

Posts: 7501
Joined: Dec 2001

Thu Aug 12, 10 03:09 PM
User is offline View users profile

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.

Edited: Thu Aug 12, 10 at 04:44 PM by Alan
 
Reply
   
Quote
   
Top
   
Bottom
     



vandervolt
Member

Posts: 51
Joined: Apr 2007

Fri Aug 13, 10 04:08 AM
User is offline

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).
 
Reply
   
Quote
   
Top
   
Bottom
     



MJPL
Junior Member

Posts: 7
Joined: Aug 2010

Fri Aug 13, 10 01:32 PM
User is offline

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
 
Reply
   
Quote
   
Top
   
Bottom
     



Cuchulainn
Senior Member

Posts: 38011
Joined: Jul 2004

Fri Aug 13, 10 03:58 PM
User is offline

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?

Thanks in advance

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
 
Reply
   
Quote
   
Top
   
Bottom
     



MJPL
Junior Member

Posts: 7
Joined: Aug 2010

Fri Aug 13, 10 04:35 PM
User is offline

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?
 
Reply
   
Quote
   
Top
   
Bottom
     



Cuchulainn
Senior Member

Posts: 38011
Joined: Jul 2004

Sat Aug 14, 10 04:55 PM
User is offline

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
 
Reply
   
Quote
   
Top
   
Bottom
     



MJPL
Junior Member

Posts: 7
Joined: Aug 2010

Sat Aug 14, 10 05:28 PM
User is offline

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
 
Reply
   
Quote
   
Top
   
Bottom
     



Cuchulainn
Senior Member

Posts: 38011
Joined: Jul 2004

Tue Aug 17, 10 01:41 PM
User is offline

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
 
Reply
   
Quote
   
Top
   
Bottom
     



Alan
Senior Member

Posts: 7501
Joined: Dec 2001

Tue Aug 17, 10 02:02 PM
User is offline View users profile

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.


Edited: Tue Aug 17, 10 at 02:07 PM by Alan
 
Reply
   
Quote
   
Top
   
Bottom
     



MJPL
Junior Member

Posts: 7
Joined: Aug 2010

Tue Aug 17, 10 02:16 PM
User is offline

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
 
Reply
   
Quote
   
Top
   
Bottom
     



Alan
Senior Member

Posts: 7501
Joined: Dec 2001

Tue Aug 17, 10 03:55 PM
User is offline View users profile

try y = x/(x+p)
 
Reply
   
Quote
   
Top
   
Bottom
     



Cuchulainn
Senior Member

Posts: 38011
Joined: Jul 2004

Wed Aug 18, 10 01:34 PM
User is offline

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
 
Reply
   
Quote
   
Top
   
Bottom
     

Pages: [ 1 2 >> Next ]
View thread in raw text format
FORUMS > Numerical Methods Forum < refresh >

Forum Navigation:

© All material, including contents and design, copyright Wilmott Electronic Media Limited - FuseTalk 4.01 © 1999-2014 FuseTalk Inc. Terms & Conditions