Diffusion Equation

Let's first define our diffusion solution.  I'll use a lower case "d" rather than an upper case "D", because the latter is the function for differentiation in Mathematica.

f=(1/(2 (Sqrt[π d t]))) Exp[((- (x^2))/(4 d t))];

Let's also define a version where we  set d*t  = τ, so that we don't have to worry about constants.

fscaled=(1/(2 (Sqrt[π τ ]))) Exp[((- (x^2))/(4 τ))];

(a)  Plot the function at several different times on the same graph.

There are a few ways to do this in Mathematica.  One is to substitute in the scaled function above. Mathematica has a way to do symbolic substitution:

f1plot=fscaled/.τ->0.1`;
f2plot=fscaled/.τ->0.5`;
f3plot=fscaled/.τ->1;
f4plot=fscaled/.τ->2;
Plot[{f1plot,f2plot,f3plot,f4plot},{x,-3,3}]

A second way is to define what we would normally consider a "function" in that we can give it an argument.  In Mathematica this is called "delayed assignment" since the function is not evaluated until you actually give it numbers.   The format is:

ffxn[x_,τ_] := (1/(2 (Sqrt[π τ ]))) Exp[((- (x^2))/(4 τ))];

The underscores and the ":=" are all part of creating the function.  Now you can invoke this as you might normally evaluate a function:

ffxn[1,3]

We can also now plot it for different values of τ.

Plot[{ffxn[x,0.1`],ffxn[x,0.5`],ffxn[x,1]},{x,-3,3}]

(b)  Verify that this is a solution of the diffusion equation.

Let's use our original f.  Below we evaluate the two derivatives, and  you can see that that they are the same up to an overall constant.  The third line just verifies the diffusion equation:

         d(( (∂^2) f)/(∂(x^2)))= (( ∂ f)/(∂t ))

D[f,x,x]//Simplify
D[f,t]//Simplify
d  D[f,x,x]== D[f,t]//Simplify

(c) Verify that the equation is normalized.

Integrate[f,{x,-Infinity,Infinity}]

This is a mess because Mathematica doesn't know if the variables d and t are positive.   Let's tell them that they are and ask again.

Integrate[f,{x,-Infinity,Infinity},Assumptions -> Re[(1/(d t))]>0 ]

Mathematica is still uncertain how to simplify this.  If we assure it that all roots can be treated as positive numbers then it will be happy:

PowerExpand[Integrate[f,{x,-Infinity,Infinity},Assumptions -> Re[(1/(d t))]>0 ] ]

(d) Calculate the average value of x and (x^2).   

Now that we have a version that works, let's just re-use it.

averageX= PowerExpand[Integrate[ f*x ,{x,-Infinity,Infinity},Assumptions -> Re[(1/(d t))]>0 ] ]
averageX2= PowerExpand[Integrate[ f*(x^2) ,{x,-Infinity,Infinity},Assumptions -> Re[(1/(d t))]>0 ] ]

(e)   We know from the handout that

           J(x,t)  = - d((  ∂f(x,t)  )/(∂x))

jfxn = -d   D[f,x]

Again we can make Mathematica happier by telling it that all roots are positive:

jfxn = PowerExpand[-d   D[f,x]]

Without loss of generality, let's set d = 1 for our plots.

j1plot= jfxn/. d-> 1 /. t-> 1

Here's the plot.  Note that the current is positive to the right of the maximum (material diffuses to the right)  and negative to the left of the maximum  (material diffuses to the left).

Plot[j1plot,{x,-3,3}]