I am new to Mathematica and I'm trying to plot the following analytical function at few time steps (say t=0.1, 0.5, 1, 10) all in one plot (similar to hold on
in MATLAB).
$$\begin{align*} p(x;t)=\frac1{\sigma\sqrt{2\pi t}}&\left[\exp\left(-\frac{(x-x_0-\mu t)^2}{2\sigma^2 t}\right)+\exp\left(\frac{-4x_0\mu t-(x+x_0-\mu t)^2}{2\sigma^2 t}\right)+\right.\\&\left.\frac{2\mu}{\sigma^2}\exp\left(\frac{2\mu x}{\sigma^2}\right)\left\{1-\Phi\left(\frac{x+x_0+\mu t}{\sigma \sqrt t}\right)\right\}\right] \end{align*}$$
This function is from page 224 (with some typos corrected) of the book:
Cox, D. R., and H. D. Miller. The Theory of Stochastic Processes. Vol. 134. CRC Press, 1977.
Note that the function $\Phi[(x+x_0+\mu t)/(\sigma \sqrt t)]$ is the integral from $-\infty$ to $x$ as:
$$\Phi(x)=\frac1{\sqrt{2\pi}}\int_{-\infty}^x e^{-\frac12 y^2}\mathrm dy$$
I am aware of Mathematica's error function Erf
and Normal Distribution but they integrate from $0$ to $x$ and my function is from $-\infty$ to $x$.
Here is what I have for a code as an initial one-time trial with some parameters.
ClearAll[μ, σ, t, z, z0, f1, f2, f3, f4, f5, f6, f];
z0 = 0.01; μ = 0.1; σ = 1; t = 10;
f1 := 1/(σ*Sqrt[2*π*t])
f2[z_] := Exp[-((z - z0 - μ*t)^2)/(2*t*σ^2)];
f3[z_] := Exp[(-4*z0*μ*t - (z + z0 - μ*t)^2)/(2*t*σ^2)];
f4[z_] := (2*μ)/σ^2*Exp[(2*μ*z)/σ^2];
f5[z] := Integrate[Exp[-1/2 ((u + z0 + μ*t)/(σ Sqrt[t]))^2]/Sqrt[2 π], {u, -∞, z}]
(*f6[z_]:=Simplify[f5,t>0];*)
f = f1[z]*(f2[z] + f3[z] + f4[z]*(1 - f5[z]))
(1/(10 Sqrt[2 π]))[z] (E^(1/200 (-0.4 - (-9.99 + z)^2)) + E^(- (1/200) (-10.01 + z)^2) + 0.2 E^(0.2 z) (1 - ((5. - 1.32904*10^-15 I) + (0.353553 (10.01 + 1. z) Erf[0.0707107 Sqrt[(10.01 + 1. z)^2]]) / Sqrt[(0.707814 + 0.0707107 z)^2])[z]))
Plot[{f[z]}, {z, z0, 10}]
which is not working (I use z
for x
here). The function should be real, but my answer has an imaginary part.
Erf[]
andNormalDistibution[]
are built-in, which may help simplify your code. $\endgroup$Erf[]
is supported:Φ[x_] := Erf[-∞, x/Sqrt[2]]/2
$\endgroup$