I am trying to use a finite difference scheme to numerically solve sixth order parabolic equations such as

\begin{equation} u_t = u_{xxxxxx} \end{equation}

with symmetry conditions \begin{equation} u_x(0,t) = u_{xxx}(0,t) = u_{xxxxx}(0,t) = 0. \end{equation}

On the vast majority of grid points, the central difference scheme \begin{equation} u_{xxxxxx} \approx \frac{1}{h^6}\left(u(x+3h) -6u(x+2h) + 15u(x+h) - 20u(x)+15u(x-h) - 6u(x-2h)+u(x-3h)\right) \end{equation} is second order accurate ($O(h^2)$) and is perfectly fine for my purposes, however I am unsure how to deal with evaluating this approximate sixth derivative at the first three gridpoints, $u_1, u_2$ and $u_3$. Extending the domain with 'fictional' gridpoints and using them to impose the symmetry conditions, e.g. \begin{equation} u_x \approx \frac{1}{2h}(u_2 - u_0) = 0,\quad u_{xxx} \approx \frac{1}{2h^3}(u_3 -2u_2 + 2u_0 - u_{-1}) = 0, u_{xxxxx} \approx \frac{1}{2h^5}(u_4 -4u_3 + 5u_2 -5u_0 +4u_{-1} -u_{-2})=0 \end{equation} allows us to substitute in the relations $u_0 = u_2$, $u_{-1}=u_3$, $u_{-2}=u_4$ into the central difference approximation around $u_1$, giving \begin{equation} u_{1xxxxxx} \approx 2u_4 -12u_3 +30u_2 -20u_1. \end{equation}

However, trying to make the same substitutions for the finite difference of $u_2$ and $u_3$ has thrown up some truly bizarre and pathological behaviour in my simulations around the origin. I have also tried ignoring any symmetry conditions at these two points and simply using the sixth order forward difference approximations, with scarcely better results.

How to proceed? Should I try somehow mixed finite difference approximations, and if so how do I factor the symmetry conditions into them?

Thanks very much.

  • $\begingroup$ You haven't said what the initial conditions are. The symmetry imposed at $x=0$ seems equivalent to using even extension of the initial values. If the initial values are like $u(x,0)=x$, the even extension is not once differentiable, which I imagine will be a problem for a PDE that involves 6th order derivative. $\endgroup$ – user147263 Sep 23 '14 at 4:20
  • $\begingroup$ You are quite right, I am interested in even, integrable solutions but posing them on the half line is simpler computationally. The initial conditions I'm working with so far have all been rescaled Gaussians of the form $A\exp(-by^2)$, so no worries about their differentiability. $\endgroup$ – Baron Mingus Sep 23 '14 at 13:03
  • $\begingroup$ I would try to do it on the entire line, for testing purpose. Would the same weirdness appear? I don't see why it would. Then you could consider why the symmetry conditions you have fail to give the same result. $\endgroup$ – user147263 Sep 23 '14 at 13:20
  • $\begingroup$ I just tried this myself, using Scilab, with initial values $\exp(-x^2)$ on $[0,1]$. I put the symmetry conditions on both sides of the interval. This is the output -- solutions gradually tends to zero. Note that I chose very small timeStep in the equation (it had to be small compared to the 6th power of spaceStep). Maybe this is the issue with your code? $\endgroup$ – user147263 Sep 25 '14 at 1:08

Your Answer

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

Browse other questions tagged or ask your own question.