NotebookDirectory[]

NotebookFileName[]

The objective is to solve the heat equation 

((∂u)/(∂t))(x,t)=(((∂^2)u)/(∂(x^2)))(x,t)  on  {(x,t) ∈  | 0≤x ≤L, 0 ≤t},

subject to the conditions

u(0,t) = 0, u(L,t)=0,  t≥0    (these are the Dirichlet BCs) 

u(x,0) = f(x), x∈[0,L]   (this is IC)

As discussed in class the sequence of product solutions which satisfy the BCs is, with n∈ℕ,

Exp[-(((((n Pi)/L)))^□)t]Sin[((n Pi)/L)x]

Then, to satisfy the IC, we form an infinite superposition of the above sequence of solutions with the Fourier coefficients which are integrals

(∑_(k=1)^∞)((2/L)(∫_0^L)f(ξ)Sin[((k Pi)/L)ξ]dξ)Exp[-(((((k Pi)/L)))^2)t]Sin[((k Pi)/L)x]

Now I implement this solution in Mathematica using symbolic integration. For simplicity, I assume that we always work on the interval [0,π]. I cannot work with infinite series, so the summation ends at some positive integer nn.  In this way I obtain an approximation for the solution;  larger nn, better approximation.

Clear[uue];uue[x_,t_,ff_,nn_]:=Module[{fcs},fcs=Table[(2/Pi)Integrate[ff[x]*Sin[k x],{x,0,Pi}],{k,1,nn}]; Sum[fcs[[k]] Exp[-(k^2)t]Sin[k x],{k,1,nn}]]

The above formula defines a function uue  of four variables. The variables x and t are the space and the time from the partial differential equation, the variable ff is the function from the initial condition and nn is a positive integer stating how many terms in the partial sum of the series we want to use. Larger the positive integer nn better the accuracy, but slower the process of calculation.  A tricky part here is the variable ff; it is a function and it has to be written as a pure function.

Pure function is a Mathematica construct in which a function is written with a “placeholder” variable # and the expression for the function ends with &. For example the square function is written as a pure function as

(((#)^2))&

To use this function we enclose the independent variable in square brackets as we always do with functions in Mathematica.

(((#)^2))&[3]

(((#)^2))&[a]

Pay attention how Pure Function is used in examples. For example for the constant function 1& (there is no variable # since it is a constant function) and with nn=10 we have

uue[x,t,1&,20]

We will explore this function later.

For a complicated function f(x) it might be slow or imposable to do symbolic integration, so I implemented the same formula using numeric integration.

Clear[uun];uun[x_,t_,ff_,nn_]:=Module[{fcsn},fcsn=Table[(2/Pi)NIntegrate[ff[x]*Sin[k x],{x,0,Pi},WorkingPrecision->20,AccuracyGoal->8,PrecisionGoal->12,MaxRecursion->60],{k,1,nn}];Chop[Sum[fcsn[[k]]*Exp[-(k^2)t]Sin[k x],{k,1,nn}]]
]

For example for the function (((sin(x)))^(1/5)) and nn=10 we have

uun[x,t,(((Sin[#]))^(1/5))&,10]

To explore the solutions we need to name each one of them with a special name, as follows (in the next example I name the solution uu0 and I keep the variables x and t since this is an approximation for the solution for a specific f(x). I do this naming with = sign since using the function uue with four variables would be very slow.

uu0[x_,t_]=uue[x,t,1&,20]

Now we explore visually how well the Fourier series approximates the function 1

Plot[{1,uu0[x,0]},{x,0,Pi},PlotRange->{0,1.3},PlotStyle->{{Thickness[0.01],Black},{Thickness[0.005],RGBColor[1,0,0]}}]

Plot3D[{uu0[x,t]},{x,0,Pi},{t,0.,2},PlotPoints->30,PlotRange->{0,1.3}]

Manipulate[Plot[{N[uu0[x,t]]},{x,0,Pi},PlotPoints->30,PlotRange->{0,1.3}],{{t,0.0125},0,3},ControlPlacement->Top]

(#(Pi-#))&[2]

Let us explore a function that satisfies the boundary conditions, like x(Pi-x)

uu1[x_,t_]=uue[x,t,#(Pi-#)&,20]

Visually explore how good is the approximation (The given function is the black thicker graph, the approximation is the thinner yellow graph; the approximation is so visually so good that it is difficult to distinguish the graphs.)

Plot[{x(Pi-x),uu1[x,0]},{x,0,Pi},PlotRange->{0,3},PlotStyle->{{Thickness[0.01],Black},{Thickness[0.004],RGBColor[1,1,0]}}]

However, when I plot the difference of the given function and its approximation we can see better what the error is.

Plot[x(Pi-x)-uu1[x,0],{x,0,Pi},PlotRange->Automatic,PlotStyle->{{Thickness[0.007],RGBColor[1,0,0]}}]

Plot3D[{uu1[x,t]},{x,0,Pi},{t,0,5},PlotRange->{0,3}]

Now explore a more complicated function which satisfies the boundary conditions:

uu3[x_,t_]=uun[x,t,(((Sin[#]))^(1/5))&,50]

Plot[{(((Sin[#]))^(1/5))&[x],uu3[x,0]},{x,0,Pi},PlotRange->{0,1.2},PlotStyle->{{Thickness[0.01],Black},{Thickness[0.004],RGBColor[1,0,0]}}]

Plot3D[{uu3[x,t]},{x,0,Pi},{t,0,2},PlotRange->{0,1.2}]