# Accuracy of numerical solutions provided by pdetoode

I am interested in solving the following differential equation using pdetoode (you can find pdetoode here).

ClearAll[PowerSubdivide]
PowerSubdivide[stop_, nfin_] := PowerSubdivide[1, stop, nfin]
PowerSubdivide[start_, stop_, nfin_] :=
Exp[Subdivide[Log[start], Log[stop], nfin]]
St = 100;
T0[x_] := Exp[-(x^2/Pi)] - x Erfc[x/Sqrt[Pi]];
domain = {lb, rb} = {0, 35};
tend = 7;
With[{T = T[x, t]}, {eq, ic,
bc} = {D[T, t] == D[T, {x, 2}] + St (1 + term) D[T, x],
T == T0[x] /. t -> 0, {T == 1 /. x -> lb, T == 0 /. x -> rb}};
termvalue = D[T, x] /. x -> lb];

points = 1000; difforder = 4; grid = Array[# &, points, domain];
ptoofunc = pdetoode[T[x, t], t, grid, difforder];
del = #[[2 ;; -2]] &;

ode = del@ptoofunc@eq /. term -> ptoofunc@termvalue;
odeic = ptoofunc@ic // del;
odebc = ptoofunc@bc;

sollst = NDSolveValue[{ode, odeic, odebc},
T /@ grid, {t, 0, tend}];(*//AbsoluteTiming*)
sol = rebuild[sollst, grid, 2];


I would like to vary the parameter St from $$0.01$$ to $$100$$. However I noticed a bad convergence of the algorithm for large St, as for St=100, expecially at low values of t. I tried to check this convergenze calculating

Derivative[1, 0][sol][0, 0](*should be -1*)
St (1 + Derivative[1, 0][sol][0, 0])(*should be 0*)


but I got unsatisfactory result for the second value for St=100, which should be $$0$$.

Furthermmore St (1 + Derivative[1, 0][sol][0, t]) does not tends to its asymptode for t -> 0 as shown in this figure:

numb1 = 20;
data1 = Table[{t , St (1 + Derivative[1, 0][sol][0, t])}, {t,
PowerSubdivide[0.000001, tend, numb1]}]
Show[{ListPlot[data1, ScalingFunctions -> {"Log", "Log"},
PlotRange -> All],
Plot[2/Pi St ArcTan[Sqrt[(4 t)/Pi]], {t, 0.000001, tend},
PlotStyle -> Red, ScalingFunctions -> {"Log", "Log"}],
Plot[St/(1 + St), {t, 0.1, tend}, PlotStyle -> {Red},
ScalingFunctions -> {"Log", "Log"}]}, Frame -> True,
PlotRange -> All]


Increasing the points does not seem to give better results. On the contrary for St=0.3, for example, the t->0 asymptode is reached, even if the data detach from it at very low values of t; for St=0.3 use

domain = {lb, rb} = {0, 45};
tend = 110;


For example for St=0.3, using 50, 100, 500, 1000 and 2000 points, I got for

St (1 + Derivative[1, 0][sol][0, 0])


-0.0105737

-0.0010041

-4.54256*10^-7

-1.4297*10^-8

-4.4702*10^-10

and

As you see, the calculated data are closer to the asymptode for t->0

2/Pi St ArcTan[Sqrt[(4 t)/Pi]]


for smaller times with the increase of the points, even if at very small times they start to detach from the asymptode due to some numerical instability.

Update After the suggestions of Alex Trounev, I tray to clarify my point. Using the algorithm proposed by Alex Trounev for St=10, but also using pdetoode, in this case rb=45 is a good approximation for Infinity:

St = 10;
T0[x_] := Exp[-(x^2/Pi)] - x Erfc[x/Sqrt[Pi]];

With[{T = T[x, t]}, {eq, ic,
bc} = {D[T, t] == D[T, {x, 2}] + St (1 + term) D[T, x],
T == T0[x] /. t -> 0, {T == 1 /. x -> lb, T == 0 /. x -> rb}};
termvalue = D[T, x] /. x -> lb];

domain = {lb, rb} = {0, 45};
tend = 10; points = 500;
difforder = 2; ugrid = Array[# &, points, domain];

m2 = NDSolveFiniteDifferenceDerivative[Derivative[2], ugrid,
DifferenceOrder -> difforder]["DifferentiationMatrix"]; m1 =
NDSolveFiniteDifferenceDerivative[Derivative[1], ugrid,
DifferenceOrder -> difforder]["DifferentiationMatrix"]; var =
Table[u[i], {i, Length[ugrid]}]; vart =
Table[u[i][t], {i, Length[ugrid]}]; dvart =
Table[u[i]'[t], {i, Length[ugrid]}]; uxx = m2.vart; ux = m1.vart;

eqs = Table[
dvart[[i]] == uxx[[i]] + St (1 + ux[[1]]) ux[[i]], {i, 2,
Length[ugrid] - 1}]; eq12 = {dvart[[1]] == 0,
dvart[[Length[ugrid]]] == 0}; ic1 =
Table[u[i][0] == T0[ugrid[[i]]], {i, Length[ugrid]}];

sol = NDSolve[{Join[eqs, eq12], ic1}, var, {t, 0, tend}];


if I plot

 Show[{Plot[St (1 + ux[[1]] /. sol[[1]]), {t, 10^-14, tend},
ScalingFunctions -> {"Log", "Log"}, PlotRange -> All],
Plot[2/Pi St ArcTan[Sqrt[(4 t)/Pi]], {t, 10^-14, tend},
PlotStyle -> {Dashed, Red}, ScalingFunctions -> {"Log", "Log"}],
Plot[St/(1 + St), {t, 0.001, tend}, PlotStyle -> {Red, Dashed},
ScalingFunctions -> {"Log", "Log"}]}, Frame -> True,
PlotRange -> All, FrameLabel -> {"t", "St (1+ux[[1]]/.sol[[1]])"},
PlotLabel -> "points=500"]


I got different results varying points. In particular, you see that increasing the points I finally reach and stay on the asymptode (see please the knee in the green circle). But after that, the solution departs from the asymptode at lower values of t. For St=100, I need points=20000 to reach the asymptode but soon after the solution detaches from the asymptode and I do not manage to stay on it. I suspect that the departure from the asymptode at lower values of t is due to numerical errors. Is it possible to correct somehow the algorithm so that the solution remains on the asymptode? Maybe without using so many points? Or maybe we need to increase the number of grid points along t, while it seems that points in Trounev's algorithm control just the the number of grid points in x?

Any help is really welcome.

• Please be more specific. 1. St (1 + Derivative[1, 0][sol][0, 0])(*should be 0*) evaluates to -0.00043487, which doesn't seem to be that bad. What accuracy are you expecting? 2. "St (1 + Derivative[1, 0][sol][0, t]) does not tends to its asymptode for t -> 0" What's the asymptode? And why do you believe the expression should tend to that asymptode? Commented Oct 6, 2023 at 14:50
• Anyway, try SetOptions[ListInterpolation, InterpolationOrder -> difforder] before executing that rebuild[...] line. Commented Oct 6, 2023 at 14:53
• No substantial change with SetOptions.
– umby
Commented Oct 6, 2023 at 19:20
• With SetOptions, St (1 + Derivative[1, 0][sol][0, 0])(*should be 0*) becomes -1.3586*10^-6 for your setting. And, please answer the questions in my first comment. Commented Oct 7, 2023 at 2:14
• Is the parameter St=100 used in that article? And, you haven't answered the other questions in my first comment so far. Commented Oct 10, 2023 at 0:11

As we know pdetoode can be used to discretize PDE, and therefore it can be replaced with any other FDM algorithm, for example

St = 100;
T0[x_] := Exp[-(x^2/Pi)] - x Erfc[x/Sqrt[Pi]];

With[{T = T[x, t]}, {eq, ic,
bc} = {D[T, t] == D[T, {x, 2}] + St (1 + term) D[T, x],
T == T0[x] /. t -> 0, {T == 1 /. x -> lb, T == 0 /. x -> rb}};
termvalue = D[T, x] /. x -> lb];

domain = {lb, rb} = {0, 35};
tend = 7; points = 4000; difforder = 2; ugrid =
Array[# &, points, domain];

m2 = NDSolveFiniteDifferenceDerivative[Derivative[2], ugrid,
DifferenceOrder -> difforder]["DifferentiationMatrix"]; m1 =
NDSolveFiniteDifferenceDerivative[Derivative[1], ugrid,
DifferenceOrder -> difforder]["DifferentiationMatrix"]; var =
Table[u[i], {i, Length[ugrid]}]; vart =
Table[u[i][t], {i, Length[ugrid]}]; dvart =
Table[u[i]'[t], {i, Length[ugrid]}]; uxx = m2 . vart; ux = m1 . vart;

eqs = Table[
dvart[[i]] == uxx[[i]] + St (1 + ux[[1]]) ux[[i]], {i, 2,
Length[ugrid] - 1}]; eq12 = {dvart[[1]] == 0,
dvart[[Length[ugrid]]] == 0}; ic1 =
Table[u[i][0] == T0[ugrid[[i]]], {i, Length[ugrid]}];

sol = NDSolve[{Join[eqs, eq12], ic1}, var, {t, 0, tend}];


If we plot in log-log coordinate function St (1+ux[[1]]/.sol[[1]]) together with asymptotic solution, then we have

Show[{Plot[St (1 + ux[[1]] /. sol[[1]]), {t, 10^-14, tend},
ScalingFunctions -> {"Log", "Log"}, PlotRange -> All],
Plot[2/Pi St ArcTan[Sqrt[(4 t)/Pi]], {t, 10^-14, tend},
PlotStyle -> Red, ScalingFunctions -> {"Log", "Log"}],
Plot[St/(1 + St), {t, 0.1, tend}, PlotStyle -> {Red},
ScalingFunctions -> {"Log", "Log"}]}, Frame -> True,
PlotRange -> All]


It looks like a large discrepancies with asymptotic solution, but this is kind of illusion due to log-log coordinates. If we plot solution in unscaled coordinates then it looks different

Show[Plot[St (1 + ux[[1]] /. sol[[1]]), {t, 0, .01},
PlotRange -> All],
Plot[2/Pi St ArcTan[Sqrt[(4 t)/Pi]], {t, 0, .01},
PlotStyle -> {Red, Dashed}, PlotRange -> {0, St/(1 + St)}],
Plot[St/(1 + St), {t, 0, .01}, PlotStyle -> {Red, Dashed}],
Frame -> True]

Show[Plot[St (1 + ux[[1]] /. sol[[1]]), {t, 0, tend},
PlotRange -> All],
Plot[2/Pi St ArcTan[Sqrt[(4 t)/Pi]], {t, 0, .01},
PlotStyle -> {Red, Dashed}, PlotRange -> {0, St/(1 + St)}],
Plot[St/(1 + St), {t, 0, tend}, PlotStyle -> {Red, Dashed}],
Frame -> True]


ux[[1]] /. sol[[1]] /. t -> 0

(*-1.*)


and

St (1 + ux[[1]] /. sol[[1]] /. t -> 0)

(*Out[]= 6.79265*10^-6 *)


Update 1. To plot solution at a given t={0,1,2,3,4,5,6,7} we use

ListLinePlot[
Table[Transpose[{ugrid, vart /. sol[[1]]}], {t, 0, 7, 1}],
PlotLegends -> Range[0, 7], PlotRange -> All]


In a small scale it looks like

ListLinePlot[
Table[Transpose[{ugrid, vart /. sol[[1]]}], {t, 0, 7, 1}],
PlotLegends -> Range[0, 7]]


• Is it possible to plot the solution, T, as function of x for a given value of t, I have not really understood the structue of the solution, sol, it looks like an array of interpolating functions at t=0.
– umby
Commented Oct 9, 2023 at 13:01
• @umby See Update 1 to my answer. Commented Oct 9, 2023 at 14:18
• @AlexTrunev, it seems that the variable points = 4000 controls only the grid spacing along x. Is it possible to control also the grid spacing along t or NDSolve takes care automatically of it? Probably increasing the numbers of points along t will produce a lower discrepancy of St (1 + ux[[1]] /. sol[[1]] with respect to the asymptote for t->0. Many thanks, also for what you have done so far.
– umby
Commented Oct 16, 2023 at 13:03