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.

Improve this question
$\endgroup$
9
  • $\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.'s eventual burnout ♦
    Apr 7, 2018 at 21:35
  • $\begingroup$ Also, note that Erf[] and NormalDistibution[] are built-in, which may help simplify your code. $\endgroup$
    – J. M.'s eventual burnout ♦
    Apr 7, 2018 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.'s eventual burnout ♦
    Apr 7, 2018 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, 2018 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.'s eventual burnout ♦
    Apr 7, 2018 at 22:03

1 Answer 1

Reset to default
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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks Henrik. This is a Probability density function and hence should be between 0 and 1. $\endgroup$
    – Khaled
    Apr 7, 2018 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, 2018 at 23:21
  • 1
    $\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, 2018 at 23:43

Your Answer

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

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