# Is it possible to solve a differential equation with a user-defined variable mesh in NDSolve?

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.

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$$.

• 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! May 3 at 16:42
• 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. May 3 at 16:50
• 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. May 3 at 17:01
• Are you simulating a liquid film under the influence of surface tension and Van der Waal's force? May 3 at 17:24
• @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. May 4 at 3:31

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:
The tutorial above is about TensorProductGrid method. As to FiniteElement, as mentioned by user21 in this answer:
Additionally, your PDE is 4th order in $$x$$ direction, while currently FiniteElement method doesn't support PDE whose spatial derivative order exceeds 2.