1
$\begingroup$

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_]:=NDSolve`FiniteDifferenceDerivative[Derivative[1],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= Developer`ToPackedArray[M, Real];
Nt=10^5;
T=Table[it*dt,{it,0,Nt-1}];
U=ConstantArray[0.,{Nt, Nx}];
U[[1, All]]=f[X];
U= Developer`ToPackedArray[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]
Improve this question
$\endgroup$
8
  • 1
    $\begingroup$ 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. $\endgroup$
    – Daniel Huber
    Jan 14, 2022 at 15:35
  • $\begingroup$ @DanielHuber I see, thank you! Is the code running correctly for you? It is returning the two latter ListPlots empty for me $\endgroup$
    – JBuck
    Jan 14, 2022 at 16:09
  • 1
    $\begingroup$ Sorry, but I do not have the time (and nerve ): ) to wait for the loop to end. $\endgroup$
    – Daniel Huber
    Jan 14, 2022 at 16:14
  • $\begingroup$ @DanielHuber Thanks anyway! $\endgroup$
    – JBuck
    Jan 14, 2022 at 16:17
  • 1
    $\begingroup$ @JBuck Sorry I didn't check your code in detail, I only tried to find an alternative solution using "MethodOfLines" $\endgroup$
    – Ulrich Neumann
    Jan 14, 2022 at 16:47

1 Answer 1

Reset to default
5
$\begingroup$

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_] := 
  NDSolve`FiniteDifferenceDerivative[Derivative[1], 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 = Developer`ToPackedArray[M, Real];
Nt = Round[1/dt];
T = Table[it*dt, {it, 1, Nt}];
U[0] = 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]}]

Figure 1

Improve this answer
$\endgroup$
8
  • $\begingroup$ 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]]=1/(2*n^2); weights[[n]]=1/(2*n^2); weights.U[0] 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[[1]]], to compute the energy,but weights.U[0] is nowhere near the actual integral value, any clue why? $\endgroup$
    – JBuck
    Jan 19, 2022 at 18:53
  • 1
    $\begingroup$ @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][[2]], but this is also not so precise since weights1 . (Drop[U[0], -1] + Drop[U[0], 1])/2=0.0619452 , and we expecting about 0.0617043. We should increase Nx up to 2^12 to get about 0.06172. $\endgroup$
    – Alex Trounev
    Jan 19, 2022 at 23:39
  • $\begingroup$ 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[[1]]], 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? $\endgroup$
    – JBuck
    Jan 20, 2022 at 14:57
  • 1
    $\begingroup$ @JBuck Ok! It looks much better {weights . U[0], 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[0], {i, 0, Nx}]] $\endgroup$
    – Alex Trounev
    Jan 20, 2022 at 16:59
  • 1
    $\begingroup$ @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[0], {i, 0, Nt}]]. $\endgroup$
    – Alex Trounev
    Jan 23, 2022 at 2:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.