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 = NDSolve`FiniteDifferenceDerivative[Derivative[2], ugrid,
DifferenceOrder -> difforder]["DifferentiationMatrix"]; m1 =
NDSolve`FiniteDifferenceDerivative[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.
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 fort -> 0" What's the asymptode? And why do you believe the expression should tend to that asymptode? $\endgroup$SetOptions[ListInterpolation, InterpolationOrder -> difforder]before executing thatrebuild[...]line. $\endgroup$SetOptions. $\endgroup$SetOptions,St (1 + Derivative[1, 0][sol][0, 0])(*should be 0*)becomes-1.3586*10^-6for your setting. And, please answer the questions in my first comment. $\endgroup$St=100used in that article? And, you haven't answered the other questions in my first comment so far. $\endgroup$