Solve a stochastic equation analytically

Question

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!

Answer 1

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!

Later edit: To address your question after you have edited it.

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) + 
 Plus @@ (MapThread[
     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!

Answer 2

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.