# Explicit integration of PDE running very slowly

I would like to solve the PDE $$\frac{\partial u}{\partial t}=(1-x)\frac{\partial u}{\partial x}$$, in the interval $$[0,1]$$, with initial condition $$u(0,x)=e^{-218(\frac{-2}{3}-ln(1-x))^2}$$, by manually implementing a pseudospectral method with Chebyshev nodes for spatial discretization and a Taylor method for time integration. The first ListPlot appears in a matter of seconds, but after that the code keeps running forever, to the point that I do not even know if it produces the correct results at the end. Is the code correct, and if so, can it become quicker? Please let me know for any clarifications.

a=218;
b=-2/3;
x0=0;
x1=1;
u[x_,t_,a1_,b1_]:=Exp[-a1*(b1+t-Log[1-x])^2];
f[x_]:=u[x,0,a,b];
D[t,u]-(1-x)*D[x,u]

der1[x0_]:=NDSolveFiniteDifferenceDerivative[Derivative,x0,
"DifferenceOrder"->"Pseudospectral"]@"DifferentiationMatrix";
Nx=2^8;
II=IdentityMatrix[Nx];
n=Nx-1;
theta=N[Table[i*Pi/n,{i,0,n}]];
X=N[Chop[((x0+x1)/2)+((x0-x1)/2)*Cos[theta]]];
dt=N[0.01*(x1-x0)/n];
D1=der1[X];
L=(1-x)*D1;
L2=L.L;
L3=L2.L;
M = (II+dt*L+((dt^2)/2)*L2+((dt^3)/6)*L3);
M= DeveloperToPackedArray[M, Real];
Nt=10^5;
T=Table[it*dt,{it,0,Nt-1}];
U=ConstantArray[0.,{Nt, Nx}];
U[[1, All]]=f[X];
U= DeveloperToPackedArray[U, Real];
ListPlot[Transpose[{X, U[[1, All]]}], PlotRange->All, Joined->True]
Do[U[[it+1,All]]=M.U [[it,All]],{it,1,Nt-1}]
ListPlot[Transpose[{X,U[[Nt, All]]}], PlotRange->All, Joined->True]
ListPlot[Transpose[{X,U[[Nt, All]]-f[X-T[[Nt]]]}], PlotRange->All, Joined->True]

• The time consuming step is: Do[U[[it+1,All]]=M.U [[it,All]],{it,1,Nt-1}]You could try Compilâtion or ParallelEvaluate. Jan 14, 2022 at 15:35
• @DanielHuber I see, thank you! Is the code running correctly for you? It is returning the two latter ListPlots empty for me Jan 14, 2022 at 16:09
• Sorry, but I do not have the time (and nerve ): ) to wait for the loop to end. Jan 14, 2022 at 16:14
• @DanielHuber Thanks anyway! Jan 14, 2022 at 16:17
• @JBuck Sorry I didn't check your code in detail, I only tried to find an alternative solution using "MethodOfLines" Jan 14, 2022 at 16:47

There are typos in the code that x in L = (1 - x)*D1 is not defined and also in the initial data last point is Indeterminate. General structure of the code should be revisited as follows

Clear["Global*"]

Needs["Developer"]

b = -2/3; a = 218;
x0 = 0;
x1 = 1;
u[x_, t_, a1_, b1_] := Exp[-a1*(b1 + t - Log[1 - x])^2];
f[x_] := u[x, 0, a, b];
D[u, t] - (1 - x)*D[u, x]

der1[x0_] :=
NDSolveFiniteDifferenceDerivative[Derivative, x0,
"DifferenceOrder" -> "Pseudospectral"]@"DifferentiationMatrix";
Nx = 2^8;
II = IdentityMatrix[Nx];
n = Nx - 1;
theta = N[Table[i*Pi/n, {i, 0, n}]];
X = N[Chop[((x0 + x1)/2) + ((x0 - x1)/2)*Cos[theta]]];
dt = N[0.01*(x1 - x0)/n];
D1 = der1[X];
L = (II - DiagonalMatrix[X]) . D1;
L2 = L . L;
L3 = L2 . L;
M = (II + dt*L + ((dt^2)/2)*L2 + ((dt^3)/6)*L3);
M = DeveloperToPackedArray[M, Real];
Nt = Round[1/dt];
T = Table[it*dt, {it, 1, Nt}];
U = Chop[f[X] /. Indeterminate -> 0] // Quiet;


Iteration on time takes about 0.25 s

Do[U[it] = M . U[it - 1], {it, 1, Nt}]


Visualization of function and logarithm of absolute error on several steps

    Table[{ListPlot[Transpose[{X, U[i]}], PlotRange -> All,
Joined -> True, PlotLabel -> i dt],
ListPlot[Transpose[{X, Abs[U[i] - u[X, T[[i]], a, b] // Quiet]}],
PlotRange -> All, Joined -> True, ScalingFunctions -> "Log"]}, {i,
1, Nt, Round[Nt/10]}] • Thank you so much! I next tried to plot the energy error $E(t)-E(0)$ as a function of t, where the energy is $\Int__{0}^{1} (1-x) (\partial_x u)^2 \,dx \$. I used the Clenshaw-Curtis integration rule h[x_]:=(1/n)*(1-Sum[2*Cos[2*i*x]/(4*i^2 - 1), {i, (n-1)/2}]) weights=h[X]; weights[]=1/(2*n^2); weights[[n]]=1/(2*n^2); weights.U NIntegrate[f[x],{x,x0,x1}] energy=ConstantArray[0.,{Nt}]; Do[energy[[it]]=weights.(D1.U[[it,All]])^2,{it,1,Nt}] ListPlot[energy-energy[]], to compute the energy,but weights.U is nowhere near the actual integral value, any clue why? Jan 19, 2022 at 18:53
• @JBuck Please, check that weights . ConstantArray[1, Length[weights]]=0.704224, and therefore weights is not defined as the weights of the Clenshaw-Curtis rule. We can use weights1 = NIntegrateClenshawCurtisRuleData[Nx/2, 20][], but this is also not so precise since weights1 . (Drop[U, -1] + Drop[U, 1])/2=0.0619452 , and we expecting about 0.0617043. We should increase Nx up to 2^12 to get about 0.06172. Jan 19, 2022 at 23:39
• Nvm I found the mistake about the weights, they should be defined as h[theta] instead of h[X]. On the other hand, when I try to compute the energy energy=ConstantArray[0.,{Nt}]; Do[energy[[it]]=weights.(D1.U[[it,All]])^2,{it,1,Nt}] ListPlot[energy-energy[]], I get errors like Part::partd: Part specification U[[1,All]] is longer than depth of object, do you have any idea why this is happening? Jan 20, 2022 at 14:57
• @JBuck Ok! It looks much better {weights . U, NIntegrate[f[x], {x, x0, x1}]} is about {0.0617043, 0.0617043}. To calculate energy we can use Do[energy[it] = weights . (D1 . U[it])^2, {it, 0, Nt}]; ListPlot[Table[energy[i] - energy, {i, 0, Nx}]] Jan 20, 2022 at 16:59
• @JBuck Don't use energy=ConstantArray[0.,{Nt}];! Run code first to compute U and weights. Also in the last line use Do[energy[it] = weights . (D1 . U[it])^2, {it, 0, Nt}]; ListPlot[[Table[energy[i] - energy, {i, 0, Nt}]]`. Jan 23, 2022 at 2:47