For some differential equations, its solution may evolve to a cusp, for example, singular behavior. One may want to introduce additional mesh points near the cusp to accurately follow the solution throughout the formation of the local singularity.

General question: Can we define a variable mesh scheme in NDSolve to solve a nonlinear PDE about $h(x,t)$ with finite difference method (FDM), in which the number of mesh points near a cusp may be determined by the value of the local minimum or maximum, $h_{min}$ or $h_{max}$. Let us consider the case of a cusp having $h_{min}$. I noted that there are three additional processes should be added: (1) after each time step, a new set of mesh needs to be defined with additional points introduced; (2) the values of $h(x^{n+1}_j,t^{n+1})$ on the new mesh needs to be obtained from those $h(x^n_i,t^n)$ on the old mesh by interpolating; (3) the time step should also be adjusted according to increasing mesh resolution to ensure a consistent evolution in local mesh scales.

For example, consider the following PDE

$$h_t+(h^3h_{xxx}+h^{-1}h_x)_x=0,$$ which is defined on $x\in[0, L]$ with $L=8$, subject to the boundary conditions $h_x=h_{xxx}=0$ at the two end points of the interval. This equation will evolve a cusp at $x=0$ with a harmonic initial condition $h(x,0)=1-0.1\cos(2\pi x/L)$.

Pleas note that I can solve this equation using NDSolve with Method of Line and some common settings, which can be found easily on the site. However, to have more control on the numerical scheme and to obtain a more accurate solution especially near the singular point, I want to use the above mentioned variable mesh scheme, which begins with, say, $20$ points and increases to a much more dense mesh by $h_{min}=10^{-5}$. Specifically, is it possible to define a mesh in NDSolve with the following transformation $x(s)=h^2_{min}(t)\sinh(s)$, where $s(x,t)$ is the coordinate used in the numerical scheme. The mesh points are uniform in $s$ and thus the mesh spacing in $x$ is nearly uniform near $x=0$ but its spacing increases away from $x=0$. Please refer to the $\sinh$ curve.

enter image description here

After each time step, a new mesh is produced using the current value of $h_{min}(t^n)$ with additional mesh points introduced to keep the mesh spacing in $s$ roughly constant. The values of $h(x^{n+1}_j,t^{n+1})$ on the new mesh are obtained from those $h(x^n_i,t^n)$ on the old mesh by cubic interpolation, say. Also, the time step is chosen by

$$\Delta t=\Delta x_0^4/h^3_{min}\approx h_{min}^5 \Delta s^4,$$ where $\Delta x_0$ is the mesh spacing near $x=0$.

  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$
    – Michael E2
    May 3 at 16:42
  • $\begingroup$ I think NDSolve can produce a solution only in the form of an InterpolatingFunction, and InterpolatingFunction can have as a domain, a tensor-product grid, a finite element mesh, or a product of the two (the product may be constrained to a one-dimensional grid for "time" by an element mesh for "space"). So while one could program the idea you propose, I don't think you can use NDSolve by itself to accomplish it. $\endgroup$
    – Michael E2
    May 3 at 16:50
  • $\begingroup$ OTOH, you might be able to have NDSolve do the time-stepping using the plug-in framework. It would probably take some work, which I don't have time to do (sorry). But a time-integration method can store and alter data (for instance, a mesh) at each step, which is basically your plan. You could Sow[] the solution vectors at each step and recover them with Reap[]. You would have to pass a solution "StepOutput" back to NDSolve. Not sure what to suggest for that. $\endgroup$
    – Michael E2
    May 3 at 17:01
  • $\begingroup$ Are you simulating a liquid film under the influence of surface tension and Van der Waal's force? $\endgroup$
    – dearN
    May 3 at 17:24
  • $\begingroup$ @dearN correct, this equation can describe such a liquid film and can be solved easily with NDSolve, as you know. Well, could you give some suggestions on how to get a more accurate solution near the singular behavior? Thanks a lot. $\endgroup$
    – user95273
    May 4 at 3:31

1 Answer 1


This is hard (if not impossible) to achieve inside NDSolve AFAIK. I sincerely hope I'm wrong. As mentioned in the tutorial The Numerical Method of Lines, particularly in Spatial Error Estimates section:

Ideally, the spatial error estimates would be monitored over time and the spatial mesh updated according to the evolution of the solution. The problem of grid adaptivity is difficult enough for a specific type of PDE and certainly has not been solved in any general way. Furthermore, techniques such as local refinement can be problematic with the method of lines since changing the number of mesh points requires a complete restart of the ODE methods. There are moving mesh techniques that appear promising for this approach, but at this point, NDSolve uses a static grid.

The tutorial above is about TensorProductGrid method. As to FiniteElement, as mentioned by user21 in this answer:

No, adaptive mesh refinement is not implemented in version 11.3 for PDE solving with the FEM. NIntegrate with the FEM uses an adaptive mesh refinement.

Additionally, your PDE is 4th order in $x$ direction, while currently FiniteElement method doesn't support PDE whose spatial derivative order exceeds 2.

  • $\begingroup$ thank you very much for the information! So it is hopeless to solve this equation using NDSolve... Since NDSolve uses a static grid, if I use "MinPoints" -> 1000, "MaxPoints" -> 2000 in the method option, will it choose a fixed value between 1000 and 2000 as actual number of points according to a prior error estimate? $\endgroup$
    – user95273
    May 4 at 3:44
  • $\begingroup$ @user95273 Yes, exactly. $\endgroup$
    – xzczd
    May 4 at 3:52

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.