# Solve a stochastic equation analytically

I have a function $\vec{F}_i(t)$, which is unknown, but I do know it's mean $\langle \vec{F}_i(t) \rangle = \vec{0}$ and it's variance $\langle \vec{F}_i(t) \cdot \vec{F}_j(t') \rangle = 2 k_B T \gamma \delta_{ij} \delta(t-t')$.

I am having a long expression and now want to evaluate it's mean - which means eventually just that I have to replace all occurrences of $\vec{F}_i(t) \cdot \vec{F}_i(t')$ by the expression above and all other occurrences by $0$.

An example function would be: $$\vec{r}(t) \cdot \vec{F}_1(\tau) =\vec{F}_1(\tau) \cdot \left(\vec{r}_0 e^{-\frac{3 t}{\gamma }} + \int _1^t \frac{\vec{F}_1(K) e^{\frac{3 K}{\gamma }} \left(e^{-\frac{3 K}{\gamma }}+2\right)}{3 \gamma } dK \right)$$

where $\vec{r}(t)$ is the result of a DSolve command. The desired result for this simple example is $$\langle\vec{r}(t) \cdot \vec{F}_1(\tau)\rangle =\int _1^t \frac{ 2 k_B T \gamma \delta(K-\tau) e^{\frac{3 K}{\gamma }} \left(e^{-\frac{3 K}{\gamma }}+2\right)}{3 \gamma } dK$$

These functions are generated by:

  rr = (r1[t]/.First @DSolve[
{
\[Gamma] r1'[t]==(r1[t]-r2[t])+(r1[t]-r3[t]) +F1[t],
\[Gamma] r2'[t]==(r2[t]-r1[t])+(r2[t]-r3[t]) +F2[t],
\[Gamma] r3'[t]==(r3[t]-r1[t])+(r3[t]-r2[t]) +F3[t]
},{
r1,
r2,
r3
},t]);
expr1 = F1[\[Tau]] rr
expr2 = rr * rr


where I have left out some prefactors. Especially tricky are the products of two integrals in expr2 where you can get a $\delta(K[1]-K[2])$ One way to do this is of course to print the formulars and to do the replacements by hand. Another way I see is to copy the output string and parse and alter it. But there must be a Mathematica-way of doing this, right?

Thanks for any help!

• I added an example to the question above. I tried Rules and failed. Do they cover case like the one above? Functions where the name of argument is unimportant? Feb 10 '15 at 11:10
• like this F[t]*g[x]* F[K] /. F[t_] F[x_] -> DiracDelta[t - x] ? Feb 10 '15 at 12:49
• I have tried this approach, but this already failes for f1[a] Integrate[(1 + f2[b]), {b, 0, 10}] /. f1[x_] f2[y_] -> DiracDelta[x - y] Feb 10 '15 at 13:08

The first replacement is easy:

F[t]*g[x]*F[k] /. u_*F[t_] F[x_] -> DiracDelta[t - x]*u

(*  DiracDelta[k - t] g[x] *)


The second one is a bit more tricky because of the integral involved. In such a case I make all substitutions first, and integrate later:

    (f1[a] (1 + f2[b]) // Expand) /. f1[x_]*f2[y_] -> DiracDelta[x - y]

(*  DiracDelta[a - b] + f1[a]  *)


Now you can integrate. If you, indeed need to preserve the integrals signs during the operation, I can only think about something like the following:

expr = f1[a] Inactivate[Integrate[(1 + f2[b]), {b, 0, 10}], Integrate];
(ReplacePart[
expr, {{1} -> 1, {2, 1} -> expr[[1]]*expr[[2, 1]] // Expand}] /.
f1[x_]*f2[y_] -> DiracDelta[x - y])

(* Inactive[Integrate][(DiracDelta[a - b] + f1[a]), {b, 0, 10}] *)


Then you may activate the result, and it will be immediately evaluated.

Have fun!

Let us denote

expr1 = F1[\[Tau]] rr // Expand


It consists of several non-integral terms followed by several terms representing products of some factors by integrals. I do not give the answer here, since it is much too long. Just evaluate and have a look at it. Now the lists lstPos and lstPosRed:

 lstPos = Position[expr1, _Integrate]
lstPosRed = lstPos /. {a_, b_} -> a

(*  {{7, 3}, {8, 4}, {9, 3}, {10, 4}, {11, 3}, {12, 4}}  *)

(*    {7, 8, 9, 10, 11, 12}       *)


gives the positions of the integrals in the expression and the positions of terms containing the integrals. With the latter we can obtain a list of all factors in front of the integrals:

lstFact = expr1[[#]] & /@ lstPosRed /. Integrate[__, _] -> 1

(*  {F1[\[Tau]]/3, 2/3 E^((3 t)/\[Gamma]) F1[\[Tau]],
F1[\[Tau]]/3, -(1/3) E^((3 t)/\[Gamma]) F1[\[Tau]],
F1[\[Tau]]/3, -(1/3) E^((3 t)/\[Gamma]) F1[\[Tau]]}   *)


while with the former we obtain the list of all integrals:

lstInt = Map[expr1[[#[[1]], #[[2]]]] &, lstPos]


Again I give no answer because of has a too large size. But you obtain it directly by evaluation. Now, here is the rule taking care of the averaging:

rule = { A_[a_] A_[b_] -> DiracDelta[a - b], A_[a_]*B_[b_] -> 0};


After that one can easily write the operator to bring the factors under the integration sighs and further averaging:

 exprFun1=Plus @@ (MapThread[
ReplaceAll[#2, #2[[1]] -> (Expand[(#2[[1]]*#1)] /.
rule)] &, {lstFact, lstInt}] //
FullSimplify[#, {t \[Element] Reals, \[Tau] \[Element] Reals}] &)


The result is very long mainly due to the HeavisideTheta...met in it many times. In order to make it more visible I (only for the sake of easier visualization here) replace HeavisideTheta[x_]->T[x]:

    Nest[ReplaceAll[#, HeavisideTheta[x_] -> T[x]] &, exprFin1, 2]

(*  -((2 E^((3 (t - \[Tau]))/\[Gamma]) (-1 + E^((
3 \[Tau])/\[Gamma])) (-1 + 2 T[-1 + t]) T[\[Tau] - t T[1 - t] -
T[-1 + t]] T[-\[Tau] + T[1 - t] + t T[-1 + t]])/(9 \[Gamma])) + (
2 E^(-((3 \[Tau])/\[Gamma])) (-1 + E^((3 \[Tau])/\[Gamma])) (-1 +
2 T[-1 + t]) T[\[Tau] - t T[1 - t] - T[-1 + t]] T[-\[Tau] +
T[1 - t] + t T[-1 + t]])/(9 \[Gamma]) + (
2 E^((3 (t - \[Tau]))/\[Gamma]) (2 + E^((3 \[Tau])/\[Gamma])) (-1 +
2 T[-1 + t]) T[\[Tau] - t T[1 - t] - T[-1 + t]] T[-\[Tau] +
T[1 - t] + t T[-1 + t]])/(9 \[Gamma]) + (
E^(-((3 \[Tau])/\[Gamma])) (2 + E^((3 \[Tau])/\[Gamma])) (-1 +
2 T[-1 + t]) T[\[Tau] - t T[1 - t] - T[-1 + t]] T[-\[Tau] +
T[1 - t] + t T[-1 + t]])/(9 \[Gamma])      *)


However, there are still terms containing no integrals. By a direct inspection one finds that they all transform into zeros upon the averaging. However, that might be by chance, while in the general case, in principle, they might be not. Let us take care of those terms. This:

    Complement[Range[Length[expr1]], lstPosRed]

(*    {1, 2, 3, 4, 5, 6}   *)


returns the positions of terms containing no integrals. Then this makes the averaging and sum those terms up:

    Plus @@ (expr1[[#]] & /@
Complement[Range[Length[expr1]], lstPosRed] /.rule)

(*  0   *)


returning zero, as expected. Finally the whole result is given by the operator:

Plus @@ (expr1[[#]] & /@
Complement[Range[Length[expr1]], lstPosRed] /. rule) +
ReplaceAll[#2, #2[[1]] -> (Expand[(#2[[1]]*#1)] /. rule)] &, {lstFact, lstInt}] //
FullSimplify[#, {t \[Element] Reals, \[Tau] \[Element] Reals}] &)


The result is evidently in this case equal to the one already obtained above. I think something analogous you can try to do with your second expression, that is, expr2 = rr * rr.

Have fun still!

• Hey - thank you for your answer, but this is not quite what I'm looking for, I think. I have edited my question. Should the Inactivate/ReplacePart work there as well? Feb 11 '15 at 11:15

Given an expression,

s = g[x] + Integrate[f[k] h[k], k] + Integrate[f[k] p[k], k]


the contraction of it with f[t] can be obtained with

s[[Flatten[Take[Position[s, f[x_]], All, 1]]]] /. f[x_] -> c DiracDelta[t - x]
(* c*Integrate[DiracDelta[-k + t]*h[k], k] + c*Integrate[DiracDelta[-k + t]*p[k], k] *)


where c is any constant. This approach also should work with other functions than Integrate, but one can imagine more complex expressions that it cannot handle.

• Thank you for your answer - it seems to work! I extended my initial question. Do you know how to handle the case s * s, i.e. expr2 in the question? Feb 11 '15 at 11:31
• @asPlankBridge I am confident that expr2 can be handled , but I shall not have time to do so for a day or so. Sorry. I did look briefly at rr, to which a straightforward modification of my approach applies. Feb 11 '15 at 15:08
• Alright - no need for you to do the work for a random internet stranger. I will look into it. Feb 11 '15 at 16:24