How to improve the solution of the following differential equation?

Consider a PDE $$\tag 1 \frac{\partial f_{n}}{\partial t}- \frac{1}{2t}E\frac{\partial f_{n}}{\partial E} = 0$$ I choose the initial condition $$f_{n}(E, t_{0}) = e^{-E/T_{0}}$$, where $$T_{0}$$ is some constant (plus the boundary condition $$f_{n}(E_{max},t) = 0$$ for some large $$E_{max}$$).

Eq. (1) may be transformed into a trivial ODE by performing a change of variables. For the given initial condition, the $$f_{n}(E,t)$$ is $$\tag 2 f_{n}(E,t) = e^{-E/(T_{0}\sqrt{t0/t})}$$ I want to obtain this result by solving the PDE. This is my code: first, I discretize the E domain, then turn the PDE into ODE by discretizing the derivative, and then solve the set of ODEs:

imax = 500;
t0 = 0.1;
(*E grid*)
DEvalue = 0.05;
Emin = 0.01;
tmax = 4;
EStep[k_, DE_] = Emin + DE*k;
(*Discretization of the derivative*)
fnEnDerivative[i_] =
If[i != imax, (fn[i + 1][t] - fn[i][t])/DEvalue, -fn[imax][t]/
DEvalue];
Tt[t_]=0.84/Sqrt[t];
(*Table with equations, initial conditions and functions*)
InitialConditionTable =
Join[{Tg[t0] == Tt[t0]},
Table[fn[i][t0] == Exp[-(EStep[i, DEvalue]/Tt[t0])], {i, 0, imax,
1}]];
functionsTable = Table[fn[i], {i, 0, imax, 1}];
EquationsTable =
Table[fn[i]'[t] -
1/(2*t)*EStep[i, DEvalue]*
fnEnDerivative[i]== 0, {i, 0,
imax, 1}];
sol = NDSolve[{EquationsTable, InitialConditionTable},
functionsTable, {t, t0, tmax},
Method -> {"EquationSimplification" -> "Solve"}][[1]];


After obtaining the solution, I compare the behavior of the exact solution $$(2)$$ and the numerical solution. I found a discrepancy due to low precision:

fnv[t_] :=
Table[{EStep[i, DEvalue], (fn[i] /. sol)[t]}, {i, 1, imax, 1}]
BoltzmannDistrDiscrete[t_] :=
Table[{EStep[i, DEvalue], Exp[-EStep[i, DEvalue]/Tt[t]]}, {i, 1,
imax, 1}]
ListLogPlot[{fnv[0.3], BoltzmannDistrDiscrete[0.3]},
PlotRange -> {{0.003, 25}, All}, Joined -> True, Frame -> True,
ImageSize -> Large]


Improving Accuracy/Precision Goals does not help. For example, with PrecisionGoal ->30 I get this:

• Value of tmax is missing. Then, you don't directly solve the PDE with NDSolve, is it kind of practice? Anyway, 1. Rule of thumb: if something is wrong with PDE solving, aside from simple mistakes, it's almost always because of improper spatial discretization, the ODE solver should always be the last thing to tackle. 2. As to AccuracyGoal and PrecisionGoal, you may want to read this: mathematica.stackexchange.com/q/118249/1871 – xzczd Nov 6 '20 at 2:32
• @xzczd Actually he has system of integrodifferential equations shown on mathematica.stackexchange.com/questions/233726/… – Alex Trounev Nov 6 '20 at 12:04

Actually the answer depends on parameters. The best what I get just playing with difference order and rational numbers is the following code

imax = 500;
t0 = 10/11;
(*E grid*)
DE = 1/20;
Emin = 1/100;
grid = Table[Emin + DE*k, {k, 0, imax}];

fd = NDSolveFiniteDifferenceDerivative[Derivative[1], grid,
"DifferenceOrder" -> 8]; m = fd["DifferentiationMatrix"];
var = Array[f, {Length[grid]}];
vart = Table[f[i][t], {i, Length[grid]}]; fx = m.vart;

eqs = Table[
f[i]'[t] - 1/(2 t) grid[[i]] fx[[i]] == 0, {i, Length[grid]}];

Tt[t_] = 84/100/Sqrt[t];
(*Table with equations,initial conditions and functions*)
InitialConditionTable =
Table[f[i][t0] == Exp[-(grid[[i]]/Tt[t0])], {i, Length[grid]}];

tmax = 4/10; sol =
NDSolve[{eqs, InitialConditionTable}, var, {t, t0, tmax},
Method -> {"EquationSimplification" -> "Solve"}][[1]];
fnv[t_] := Table[{grid[[i]], (f[i] /. sol)[t]}, {i, Length[grid]}]
BoltzmannDistrDiscrete[t_] :=
Table[{grid[[i]], Exp[-grid[[i]]/Tt[t]]}, {i, Length[grid]}]
ListLogPlot[{fnv[tmax], BoltzmannDistrDiscrete[tmax]},
PlotRange -> {{.003, 25}, All}, Joined -> True, Frame -> True,
ImageSize -> Large]


And on the next step we can use option WorkingPrecision -> 30 to clean up plot as Bob Hanlon recommended

tmax = 4/10; sol =
NDSolve[{eqs, InitialConditionTable}, var, {t, t0, tmax},
Method -> {"EquationSimplification" -> "Solve"},
WorkingPrecision -> 30][[1]];


Finally we have desired agreement

• In the NDSolve use WorkingPrecision -> 30` to clean up the plot. – Bob Hanlon Nov 6 '20 at 7:33
• @BobHanlon Thank you very much! It is been my next step I forgot to make. – Alex Trounev Nov 6 '20 at 12:01