# Numerically solving PDE with high precision

I want to numerically solve the PDE

$\partial_t u(t,x)=c\partial_x u(t,x)+(mx-l)u(t,x)$

with some initial and boundary conditions and given parameters $c$, $m$ and $l$.

Consider the code

T=55;
xMax=200;
l=.6;
m=0.0072;
c=2;

system={
D[u[t,x],t]==c*D[u[t,x],x]+(m*x-l)*u[t,x],
u[0,x]==500*Evaluate@PDF[TruncatedDistribution[{0,xMax-1},NormalDistribution[100,10]],x],
u[t,0]==u[t,xMax]==0
};

(*prepare animation of progress update*)
"Computation started: "<>DateString[]
currentTimeStep=0;
"Time step: "ProgressIndicator[Dynamic[currentTimeStep],{0,T}]

sol=First@NDSolve[system,u[t,x],{t,0,T},{x,0,xMax},StepMonitor:>(currentTimeStep=t;),Method->{"MethodOfLines","TemporalVariable"->t,"SpatialDiscretization"->{"TensorProductGrid","MinPoints"->10000}}];

"Computation stopped: "<>DateString[]

Manipulate[
Plot[Evaluate[u[t,x]/.sol/.{t->time}],{x,0,xMax},PlotRange->{{0,xMax},{-50,50}},AxesOrigin->{0,0}]
,{time,0,T}]


As you can see, for $T=55$, the solution shows numerical artifacts.

It seems that the cause for this is a not sufficent accuracy of the computation, since

u[t,x]/.sol/.{t->10,x->#}&/@Range[xMax-5,xMax]


yields

{2.53114*10^-18, -1.23827*10^-17, -7.49676*10^-15, -2.81533*10^-15, 2.97761*10^-14, 0.}


but it follows from the PDE that the solution stays non-negative for non-negative initial data.

Using a higher AccuracyGoal, such as AccuracyGoal->6 only results in huge memory consumption, even with a broader grid-spacing.

Any ideas how to eliminate the produced artifacts?

-
I've tried setting WorkingPrecision and PrecisionGoal to 10 with 200 MinPoints and the artifacts were gone. – swish Oct 26 '12 at 19:06
You should try using the Method->{"LSODA"} unless there is a reason you use a staggered grid through TensorProductGrid. I ran your simulation and with AbsoluteTiming I got 0.06 seconds of run time with LSODA and 8+ seconds with the MethodOfLines. :) Just thought that would be interesting for you! – drN Oct 27 '12 at 2:11