SOLVING THE PDE WITH A SHARP INITIAL CONDITION

The calculation below shows how the solution to the PDE may be obtained by modifying the MaxSteps and StartingStepSize options. The graph of the initial condition in the problem
agrees well with that in the solution.

The problem

Here is an example of a differential equation where the boundary conditions suggest an abrupt change of direction for the solution function. If we try to solve this equation using the default settings of NDSolve we do not get good results.

eqd1=NDSolve[{ Derivative[0,1][P][x,t]== Derivative[2,0][P][x,t],P[x,0]==(10/(Sqrt[π]))Exp[-100x^2],P[-5,t]==P[5,t]},P,{x,-5,5},{t,0,5}]

A plot of the solution

Plot3D[P[x,y]/.eqd1[[1]],{x,-5,5},{y,0,5}]

We see that the boundary conditions are not well satisfied since the boundary function is always positive and the near edge of the above surface clearly takes negative values.

Getting a better solution

In cases such as this it is often necessary to force Mathematica to take smaller steps at the outset. We can do that by specifying the starting step size. This will usually require that we increase the maximum number of steps allowed. This does give a warning message but we will see that the solution appears to be a good one. Since we are computing many more points this will take longer to run.

eqd2=NDSolve[{ Derivative[0,1][P][x,t]== Derivative[2,0][P][x,t],P[x,0]==(10/(Sqrt[π]))Exp[-100x^2],P[-5,t]==P[5,t]},P,{x,-5,5},{t,0,5},MaxSteps->500,StartingStepSize->0.05]

A plot of the second solution

Plot3D[Evaluate[P[x,t]/.First[eqd2]],{x,-5,5},{t,0,5},PlotPoints->30,PlotRange->All]

Satisfied boundary conditions

In the following plot of the boundary function the given boundary conditions are plotted in red and the boundary obtained in the solution in black.

Plot[{(10/(Sqrt[π]))Exp[-100x^2],Evaluate[P[x,0]/.First[eqd2]]},{x,-0.5,0.5},PlotRange->All, PlotStyle -> {{AbsoluteThickness[3],Hue[.95]},GrayLevel[0]}]