I am trying to reproduce the Faynman-Kac results from this Wikipedia page: Faynman-Kac

This is the equation that I am trying to solve:

$$ \frac{\partial w(t,x)}{\partial t} = \frac{1}{2} \frac{\partial^2 w(t,x)}{\partial x^2} - uV(x)w(t,x) $$ With initial condition $ w(0,x)= \delta(x) $.

As an example lets set

$$ uV(x) = ux^2 $$

I tried two different methods. I knew there was no way the first method was not going to work - but I believe in magic. Perhaps, someone knows some trick to propagate the initial condition.

Method 1 - NDSolve

w[u_] = (-u*x^2*f[t, x] + 0.5*D[f[t, x], x, x]);

wsol[u_] := NDSolve[{
     D[f[t, x], t] == w[u],
     f[0, x] == DiracDelta[x],
     f[t, -50] == Exp[-1000 t],
     f[t, 50] == Exp[-1000 t]},
    f, {t, 0, 100}, {x, -50, 50}, MaxStepSize -> 0.5, 
    AccuracyGoal -> 8, PrecisionGoal -> 8, 
    Method -> {"MethodOfLines", 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "MinPoints" -> 1000}}];

Method 2 - Fourier Transforms

    FourierTransform[1/2*D[f[t, x], x, x], x, w]

(* -(1/2) w^2 FourierTransform[f[t, x], x, w] *)

FourierTransform[-u*x^2, x, w]

(* Sqrt[2 π] u (DiracDelta^′′)[w] *)

FourierTransform[DiracDelta[x], x, w]

(* 1/Sqrt[2 π] *)

DSolve[{D[f[t], t] == 
       f[t]*(Sqrt[2 π] u (DiracDelta^′′)[w] - 1/2 w^2 ),
       f[0] == 1/Sqrt[2 π]}, f[t], t]

(* {{f[t] -> E^(
       1/2 t (-w^2 + 2 Sqrt[2 π] u (DiracDelta^′′)[w]))/
       Sqrt[2 π]}} *)

So, at this point - I get a solution but I wish I could get the inverse.

 1/2 t (-w^2 + 
    2 Sqrt[2 π] u (DiracDelta^′′)[w]))/Sqrt[
 2 π], w, x]
(* InverseFourierTransform[E^(
 1/2 t (-w^2 + 
    2 Sqrt[2 π] u (DiracDelta^′′)[w]))/Sqrt[
 2 π], w, x] *)

I think its the transform of $ x^2 $ that Mathematica does not know how to inverse.

Now, I also tried substituting the Delta function with the derivative of the HeavisideTheta[] but that still does not work.

Can someone help please? The Faynman-Kac is a very important theorem.


Based on the comment of xzczd - I was able to get some results. However, I can't get rid of the errors when I plot. Can I get some help with the plot errors please.

Numerical method - This is with a Triangular Distribution Probability Density but I get the same errors with xzczd's function. I followed Mathematica's suggestion but they do not get resolved.

w3[u_] = (-u*x^2*f[t, x] + 0.5*D[f[t, x], x, x]);

wsol3[u_] := NDSolve[{
    D[f[t, x], t] == w3[u],
    f[0, x] == 
     Evaluate[PDF[TriangularDistribution[{-0.001, 0.001}, 0.0], x]],
    f[t, -50] == 0,
    f[t, 50] == 0},
   f, {t, 0, 100}, {x, -50, 50}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid"}}];

 Plot3D[Evaluate[f[t, x] /. wsol3[u]], {t, 0, 10}, {x, -5, 5}, 
  PlotRange -> {0, 10}], {u, 0.01, 1}]
  • 1
    $\begingroup$ NDSolve can't handle DiracDelta (at least now). A common way to circumvent this is to use a very sharp function e.g. With[{a = 5}, a/Sqrt[Pi] Exp[-(a x)^2]] instead. As far as I can tell, FourierTransform can't be used to solve this problem, because the equation contains a term with variable coefficient i.e. u x^2 w[t, x]. $\endgroup$
    – xzczd
    Jul 29, 2016 at 11:36
  • $\begingroup$ @xzczd Hey thank you. That is an awesome idea. Furthermore, the property of the delta function is that the area under the curve is one. So, in line with what you suggested, a triangular distribution would work very well - I think. You can make it extremely narrow and the area under the curve would remain one. Does it sound feasible? $\endgroup$
    – Edv Beq
    Jul 29, 2016 at 23:23
  • $\begingroup$ As to the update of the question: your initial condition and boundary condition is indeed inconsistent, the i.c. suggests f[0, -50]==0, while the b.c. suggests f[0,-50] == 1. What's the actual b.c. you're trying to set? $\endgroup$
    – xzczd
    Jul 30, 2016 at 5:18
  • $\begingroup$ @xzczd Ah yes both of the b.c. can be zero but is the other error "Using maximum number of points..." that I can't get rid of. $\endgroup$
    – Edv Beq
    Jul 30, 2016 at 7:41
  • $\begingroup$ When using a piecewisely smooth initial condition, it's not that rare to see this warning showing up, usually a little option adjusting will help, for example, {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement", "MeshOptions" -> {MaxCellMeasure -> 0.0001}}. Using the smooth i.c. I suggested can also resolve the problem. (Personally I prefer this solution. ) $\endgroup$
    – xzczd
    Jul 30, 2016 at 8:34

1 Answer 1


This problem looks interesting enough to do both analytically and numerically. Start with the analytical solution using V[x] = x^2

pde = (1/2)*D[w[x, t], x, x] - D[w[x, t], t] - u*x^2*w[x, t] == 0

Separate variables

w[x_, t_] = X[x] T[t]

pde/w[x, t] // Expand
(*-T'[t]/T[t]-u x^2 + X''[x]/(2 X[x]) == 0*)

% + T'[t]/T[t]

Each side must be equal to a constant. We will use a negative constant necessary to satisfy the bc's. Solve the x equation

xeq = X''[x]/(2 X[x]) - u x^2 == -α^2

DSolve[xeq, X[x], x] // Flatten

X1 = X[x] /. % /. {C[1] -> c1, C[2] -> c2}
(*c1 ParabolicCylinderD[(Sqrt[2] α^2 - Sqrt[u])/(2 Sqrt[u]), 
   2^(3/4) u^(1/4) x] + 
 c2 ParabolicCylinderD[(-Sqrt[2] α^2 - Sqrt[u])/(2 Sqrt[u]), 
   I 2^(3/4) u^(1/4) x]*)

teq = T'[t]/T[t] == -α^2

DSolve[teq, T[t], t] // Flatten

T1 = T[t] /. % /. C[1] -> 1
(*E^(α^2 (-t))*)

Now we have a first pass at w

w[x_, t_] = X1 T1

First off, go for real solutions by getting rid of c2.

c2 = 0

Which leaves us with:

w[x, t]
c1 E^(α^2 (-t)) ParabolicCylinderD[(Sqrt[2] α^2 - Sqrt[u])/(2 Sqrt[u]),2^(3/4) u^(1/4) x]

The problem doesn't really state bc's, but I assume we want the solution bounded at x = +-infinity which requires the first argument of the ParabolicCylinderD to be an integer, and that will set α

α = α /. Solve[(Sqrt[2] α^2 - Sqrt[u])/(2 Sqrt[u]) == n, α][[2]]
(*(Sqrt[2 n + 1] u^(1/4))/2^(1/4)*)

$Assumptions = n ∈ Integers

Now w is

w[x_, t_] = c1 E^(-α^2 t) ParabolicCylinderD[n, 2^(3/4) u^(1/4) x]

We will have an infinite sum as n goes from 0 to ∞. This solutions solves the end point bc's, so now we just need to solve c1 from the initial condition.

w[x, 0] == DiracDelta[x]
c1 ParabolicCylinderD[n, 2^(3/4) u^(1/4) x] == DiracDelta[x]

Use orthogonality to solve for c1.

c1*Integrate[ParabolicCylinderD[n, 2^(3/4)*u^(1/4)*x]^2, {x, -Infinity, Infinity}] == 
  Integrate[DiracDelta[x]*ParabolicCylinderD[n, 2^(3/4)*u^(1/4)*x], {x, -Infinity, Infinity}]

MMa cannot do the first integral for general n, but we can infer it from:

Table[Integrate[ParabolicCylinderD[n, 2^(3/4)*u^(1/4)*x]^2, {x, -Infinity, Infinity}], {n, 0, 5}]
(*{Sqrt[π]/2^(1/4), Sqrt[π]/2^(1/4), 2^(3/4) Sqrt[π], 
 3 2^(3/4) Sqrt[π], 12 2^(3/4) Sqrt[π], 
 60 2^(3/4) Sqrt[π]}*)

(c1*n!*Sqrt[Pi])/(2^(1/4)*u^(1/4)) == Integrate[DiracDelta[x]*ParabolicCylinderD[n, 2^(3/4)*u^(1/4)*x], 
   {x, -Infinity, Infinity}]

c1 = c1 /. Solve[%, c1][[1]]

w[x, t] will be an infinite sum from n = 0 to ∞, but for computing and plotting use a finite sum. Use wan for the analytic solution.

wan[x_, t_, mm_] := (2*u)^(1/4)*Sum[(2^(n/2)*ParabolicCylinderD[n, 2^(3/4)*u^(1/4)*x])/(E^(((2*n + 1)*t*Sqrt[u])/Sqrt[2])*(n!*Gamma[1/2 - n/2])), {n, 0, mm}]


The numeric solution. We obviously cannot use the actual DiracDelta for the numeric solution, but we can use an approximation that has nearly the same inpulse effect. Try

f[x_] = aa E^(-10000 x^2)

Choose aa so that like the DiracDelta, the area under the curve is 1.

Integrate[f[x], {x, -Infinity, Infinity}] == 1;

aa = aa /. Solve[%, aa][[1]]

The numeric solution. A little quicker than the analytic one. First we need to set u

u = 1

sol = NDSolve[{pde, w[x, 0] == f[x], w[100, t] == 0, w[-100, t] == 0},
    w[x, t], {x, -100, 100}, {t, 0, 10}, 
   MaxSteps -> {50000, Automatic}] // Flatten

wnum[x_, t_] = w[x, t] /. sol

We are setting x = +-100 for the boundaries instead of inf. The solution damps out much sooner than those values anyway. Comparison plots of analytic and numerical solutions

Plot[{wnum[x, .01], wan[x, .01, 200]}, {x, -2, 2}, PlotRange -> All]

enter image description here

Plot[{wnum[x, .1], wan[x, .1, 50]}, {x, -2, 2}, PlotRange -> All]

enter image description here

Both are reasonable matches for times away from 0. For smaller times the analytic solution converges very slowly and takes numerous terms. The more terms used, the higher the peak at x = 0, t = 0 but eventually with many more terms computational instability occurs. The exponential term promotes quick convergence with a higher time. Animating the numerical solution:

tb = Table[ Plot[Evaluate[wnum[x, t]], {x, -5, 5}, PlotRange -> {0, 1.5}], {t, 0, 6, .1}];

enter image description here

As requested a 3D plot

Plot3D[wnum[x, t], {x, -5, 5}, {t, 0, 5}, PlotRange -> {0, 1}, 
 AxesLabel -> Automatic]

enter image description here

  • $\begingroup$ Thank you so much for your answer. Very informational. Unfortunately, my student license expired but I would love to see a 3d plot of the solution - if you could add it to the answer it would be very fulfilling. For example from t[0, 5] and x[-5, 5] $\endgroup$
    – Edv Beq
    Dec 21, 2018 at 13:59
  • $\begingroup$ wow - this is amazing to me. I have been able to solve this, as faynman and kac suggested, by generating diffusion paths - but you can only solve for one fixed initial state condition x_0=X at a time. One can imagine how computationally expensive that calculation is. Thank you so much again. $\endgroup$
    – Edv Beq
    Dec 22, 2018 at 16:17
  • $\begingroup$ It is an interesting problem. I had never actually used ParabolicCylinder functions before. $\endgroup$
    – Bill Watts
    Dec 22, 2018 at 18:45
  • $\begingroup$ I have a similar problem @BillWatts , could you have a look if you have time? $\endgroup$ Feb 12, 2020 at 17:06
  • $\begingroup$ @rami_salazar Is the problem posted somewhere? $\endgroup$
    – Bill Watts
    Feb 12, 2020 at 23:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.