0
$\begingroup$

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}]

my attempt

which is not working (I use z for x here). The function should be real, but my answer has an imaginary part.

$\endgroup$
  • $\begingroup$ Please put your code as text that can be copied by other people (select your code in Mathematica and right-click: Copy As > Input Text), and not as pictures. Otherwise, people are less likely to help you. $\endgroup$ – J. M. is away Apr 7 '18 at 21:35
  • $\begingroup$ Also, note that Erf[] and NormalDistibution[] are built-in, which may help simplify your code. $\endgroup$ – J. M. is away Apr 7 '18 at 21:36
  • $\begingroup$ "Mathematica's error function is from $0$ to $x$." - that's why two-argument Erf[] is supported: Φ[x_] := Erf[-∞, x/Sqrt[2]]/2 $\endgroup$ – J. M. is away Apr 7 '18 at 21:49
  • $\begingroup$ Thanks for your comments. I tried copying the code as "Input Text" but it still looks weird on this screen. I wonder if it has to do with how I am copying it here as code. Thanks $\endgroup$ – Khaled Apr 7 '18 at 21:59
  • $\begingroup$ It's fine; the point is to not make people have to retype everything you did, as that decreases the chances of people helping you. $\endgroup$ – J. M. is away Apr 7 '18 at 22:03
1
$\begingroup$

I fixed some of the function definitions and changed enforced evaluation of Integrate at the time of definition with Evaluate (otherwise the Plot takes forever as the symbolic integration has to be done over and over again).

ClearAll[μ, σ, t, z, z0, f1, f2, f3, f4, f5, f6, f];
z0 = 0.01; μ = 0.1; σ = 1; t = 10;
f1[z_] := 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_] := Evaluate[Integrate[
   Exp[-1/2 ((u + z0 + μ*t)/(σ Sqrt[t]))^2]/
    Sqrt[2 π], {u, -∞, z}]]
(*f6[z_]:=Simplify[f5,t>0];*)

f[z_] := f1[z]*(f2[z] + f3[z] + f4[z]*(1 - f5[z]))

Plot[{f[z]}, {z, z0, 10}]

At least, this returns plot.

$\endgroup$
  • $\begingroup$ Thanks Henrik. This is a Probability density function and hence should be between 0 and 1. $\endgroup$ – Khaled Apr 7 '18 at 23:05
  • 1
    $\begingroup$ @Khaled Probability density functions are not limited to be between 0 and 1 (although cumulative distribution functions are limited to that range). $\endgroup$ – JimB Apr 7 '18 at 23:21
  • $\begingroup$ f5[z_] could be changed to f5[z_]:=CDF[NormalDistribution[0, 1], (z + z0 + \[Mu] t)/(\[Sigma] Sqrt[t])] or 1/2 Erfc[-((z + z0 + t \[Mu])/(Sqrt[2] Sqrt[t]))]. $\endgroup$ – JimB Apr 7 '18 at 23:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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