4
$\begingroup$

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]

enter image description here

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

Could you please help me with how to resolve this problem?

$\endgroup$
2
  • 1
    $\begingroup$ 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 $\endgroup$ – xzczd Nov 6 '20 at 2:32
  • $\begingroup$ @xzczd Actually he has system of integrodifferential equations shown on mathematica.stackexchange.com/questions/233726/… $\endgroup$ – Alex Trounev Nov 6 '20 at 12:04
6
$\begingroup$

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 = NDSolve`FiniteDifferenceDerivative[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]

Figure 1

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
Figure 2

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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