Initial condition 1

We want to visualize a solution for the heat equation on the half-line.

Start by defining the fundamental solution:

g[x_,t_]=Exp[-x^2/4/k/t]/Sqrt[4Pi k t]

Now define an initial condition:

u0[x_]=Exp[-x^2]

Plot the initial condition:

Plot[u0[x],{x,0,5},PlotRange->All]

Ask Mathematica to evaluate the relevant integral:

Integrate[u0[y](g[x-y,t]-g[x+y,t]),{y,0,Infinity},Assumptions->{k>0,t>0}]

Success so let’s define u:

u[x_,t_]=%

Pick a value for k so we can plot and animate:

k=1;
Plot[{u0[x],u[x,0.001]},{x,0,10},PlotRange->All]

k=1;
Animate[Plot[{u0[x],u[x,t]},{x,0,10},PlotRange->All],{t,0,1},AnimationRunning->False]

Initial condition 2

We want to visualize a solution for the heat equation on the half-line.

Start by defining the fundamental solution:

Clear[k]

g[x_,t_]=Exp[-x^2/4/k/t]/Sqrt[4Pi k t]

Now define an initial condition:

u0[x_]=Sin[x]/x

Plot the initial condition:

Plot[u0[x],{x,0,25},PlotRange->All]

Integrate[u0[y](g[x-y,t]-g[x+y,t]),{y,0,Infinity},Assumptions->{k>0,t>0}]

Mathematica does not seem to do anything with this integral so we need a new strategy.  Let’s try some numerical integration using NIntegrate.  For this, we need a specific value of k.

k=1;

The argument for NIntegrate can have only one variable, namely the integration variable y in this case.  So, we need to supply a specific values for x and t.  Here is the choice x=1, t=1:

NIntegrate[u0[y](g[1-y,1]-g[1+y,1]),{y,0,Infinity}]

To define u[x,t] using NIntegrate,  we need to use a “delayed assignment” with the synatax :=.  With this type of assignment, the right side is not evaluated immediately.  Instead, the evaluation is done each time the assignment is called.

u[x_,t_]:=NIntegrate[u0[y](g[x-y,t]-g[x+y,t]),{y,0,Infinity}]

Check that this works as expected by evaluating u[1,1] and comparing to the result above:

u[1,1]

k=1;
Plot[{u0[x],u[x,0.1]},{x,0,10},PlotRange->All]

The Plot function needs to repeatedly evaluate NIntegrate which is time-consuming so it takes a long time to produce this plot.  To quantify this, we can use the Timing function:

Timing[Plot[{u0[x],u[x,0.1]},{x,0,10},PlotRange->All]]

The Timing function returns a list in which the first item is the time (in seconds) and the second is the output of the function being timed.  So, this plot takes about 8 seconds to produce.  To compare, here’s the timing for a simple plot:

ClearSystemCache[]

Timing[Plot[u0[x],{x,0,10},PlotRange->All]]

Much faster!  Since producing one frame of an animation would take about 8 seconds, we won’t bother.