This is a follow-up question to this previous post by @FLP, in which an interesting system of equations was solved with the useful pdetoode developed by @xzczd. I have tried to solve this problem with a different initial condition (IC). And I appreciated the help and effort from @xzczd and @FLP. The following code is a combination of the original post and the answer. Due to the symmetry, the differential equation for h[r,t]
is to be solved from $r=0$ to a sufficiently large value $r>>1$. Here, the left bound of $r$ is taken as a small non-zero value due to some numerical reason, and the right bound is taken as $50$.
My question:
In the original post, a system about $h[r,t]$ was solved with an IC of the exponential function IC = {h[r, t] == 20 E^-(r/20)^2 + hf} /. t -> 0;
. For some weird reason if you change the amplitude of the IC, say, IC = {h[r, t] == 10 E^-(r/10)^2 + hf} /. t -> 0;
, it would stuck at a small time (relative to a given Tend
) and return an error: step size is effectively zero; singularity or stiff system suspected.
I also tried a parabolic IC defined by If
:
IC = {h[r, t] == unitStepExpand@If[r >= lb && r < 1, 3/5*(1 - r^2) + hf, hf]} /. t -> 0
but the error of a singularity remained. Noted that in this case PiecewiseExpand
transforms If[...]
to an equivalent Piecewise[...]
; then SimplifyPWToUnitStep
transforms Piecewise[...]
to an equivalent combination of UnitStep[...]
. This is a necessary transformation since pdetoode
assumes the functions in the equations, ICs, and BCs are all Listable
.
For the piecewise IC, I guess the issue may be caused by the discontinuous derivative at $r=1$. But I cannot understand why the code fails with the IC only changed in its amplitude: h[r, t] == 10 E^-(r/10)^2 + hf
.
An even more strange issue is that if I change the side of the domain with, say rb=100
, a singularity occurs again.
The code:
The code for pdetoode
neglected here, please find in the abovementioned link. Here, Int[r,t]
is an integral function, defined as
$$Int(r,t)=\int_0^\infty K(r, r^\prime)\frac{\partial}{\partial r^\prime }\left(\frac{h_f}{h} \right)^3 \mathrm{d}r^\prime, $$
where $K(r,r^\prime)=\frac{2}{\pi}r[\mathbf{K}(r^\prime/r)-\mathbf{E}(r^\prime/r)]$ for $r^\prime<r$
or $K(r,r^\prime)=\frac{2}{\pi}r^\prime[\mathbf{K}(r/r^\prime)-\mathbf{E}(r/r^\prime)]$ for $r^\prime>r$, in which $\mathbf{K}$ and $\mathbf{E}$ are the complete elliptic integrals.
In order to avoid the singularity of $K(r,r^\prime)$ at $r=r^\prime$,
$r$ takes the values at the grid point lb+i*step
and $r^\prime$ takes the values at the staggered grid lb+(i+0.5)*step
.
(*constants*)
a = 0.001; hf = 0.01; \[Mu] = 1; \[Beta] = 1;
(*Equations*)
With[{h = h[r, t], P = P[r, t], Q = Q[r, t], J = J[r, t]},
EqnP = P == D[h, r, r] + D[h, r]/r + a^2/h^3;
EqnQ = Q == 1/3*h^3*r*D[P, r];
EqnC = D[h, t] + (\[Mu] D[Q, r])/r == -J;
EqnJ = J == (\[Beta] D[Int[r, t], r])/r;]
(*Computation Domain*)
lb = 1/100; rb = 50;
(*unitStepExpand is a function that transforms If[\[Ellipsis]] to mathematically equivalent combination of UnitStep[\[Ellipsis]], useful in the computation of the following integral*)
unitStepExpand = Simplify`PWToUnitStep[PiecewiseExpand[#]] &;
(*initial and boundary conditions*)
IC = {h[r, t] == 20 E^-(r/20)^2 + hf} /. t -> 0;
(*IC = {h[r, t] == 10 E^-(r/10)^2 + hf} /. t -> 0;*)
(*IC = {h[r, t] == unitStepExpand@If[r >= lb && r < 1, 3/5*(1 - r^2) + hf, hf]} /. t -> 0;*)
BC = {{D[h[r, t], r] == 0, D[h[r, t], r, r, r] == 0, D[h[x, t], x] == 0, D[h[x, t], x, x, x] == 0} /. {r -> lb, x -> rb}};
(*Discretization setting*)
difforder = 2;
points = 200;
grid = Array[# &, points, {lb, rb}];
hgrid = Map[h[#][t] &, grid];
hdgrid = fdd[Derivative[1], grid, hgrid];
(*the Integral Function*)
Kernel[r_, i_] := (2*r)/\[Pi]*unitStepExpand@If[r > i, r (EllipticK[i/r] - EllipticE[i/r]), i (EllipticK[r/i] - EllipticE[r/i])];
Integrand = D[(hf/h[r, t])^3, r] /.Thread[{D[h[r, t], r], h[r, t]} -> {hdgrid, hgrid}];
step = (rb - lb)/points;
Int[r_][t_] := Total@Table[step*Kernel[r, lb + (i + 0.5)*step]*0.5 (Integrand[[i]] + Integrand[[i + 1]]), {i, 1, points - 1}];
(*Discretising the rest of the system*)
ptoofunc = pdetoode[{h[r, t], P[r, t], Q[r, t], J[r, t], Int[r, t]}, t, grid, difforder];
removeredundant = #[[3 ;; -3]] &;
odeadd = Map[ptoofunc, {EqnP, EqnQ, EqnJ}, {2}];
ode = Block[{P, Q, J}, Set @@@ odeadd; ptoofunc@EqnC] // removeredundant;
odIC = ptoofunc@IC;
With[{sf = 1}, odBC = diffbc[t, sf]@BC // ptoofunc];
(*Solving the System*)
var = Map[h, grid];
showStatus[status_] := LinkWrite[$ParentLink,SetNotebookStatusLine[FrontEnd`EvaluationNotebook[],ToString[status]]];
clearStatus[] := showStatus[""];
clearStatus[]
monitor[t_] := EvaluationMonitor :> showStatus["t = " <> ToString[CForm[t]]]
Tend = 200;
sollst = NDSolveValue[{ode, odIC, odBC}, var, {t, 0, Tend}, monitor[t], Method -> {EquationSimplification -> MassMatrix}]; // AbsoluteTiming
(*postprocess*)
solfunc = rebuild[sollst, grid];
Manipulate[Plot[solfunc[t, r], {r, lb, rb}, PlotRange ->All], {t, 0, Tend}]
Some notes:
This code works well with the original IC, but when you change the amplitude of the IC, say, from $20$ to $10$, a singularity will be encountered. Increasing the point number and order of FD has little improvement.
I am wondering if a nonuniform mesh (or adaptive mesh) could be implemented in the framework of
pdetoode
such that a refined resolution is bunched towards the region where the change of slope ofh[r,t]
is large based on physical reason?As mentioned by @FLP, the original code has some additional problems as well. I'm not sure if the treatment of IC will be one of the potential issues. So I post this question here to draw attention.
Please give me some suggestions. Thank you in advance!
Update:
With the suggestion of xzczd, the code can run with
IC = {h[r, t] == 10 E^-(r/20)^2 + hf} /. t -> 0;
, which only changes in its amplitude, up to T=502 approximately, and then singularity occurs at the radius position of the blob, where the slope changes most rapidly.By changing the width of the initial blob, I found smaller the initial width earlier the singularity occurs. For example, using
IC = {h[r, t] == unitStepExpand@If[r >= lb && r < 1, 3/5*(1 - r^2) + hf, hf]} /. t -> 0
, the singularity occurs at t=0.000024 approximately. Increasingpoints
improves things a little.