Is the discretization for the PDE correct?

[member 428835]

Hi PF!

I'm trying to finite difference a FTCS PDE $$h_t+\partial_x\left[h^3(h_{xxx}+(Kf'(h)-G)h_x\right] = 0$$ where $K$ and $G$ are constants. $f'(h) = -(n(h^*/h)^n-m(h^*/h)^m)/h$. Boundary conditions are $h_x=h_{xxx}=0$ at both ends of the $x$ domain (however long you want to make it) and some initial height, call it $h_0(1+\epsilon \cos(2 \pi x/\lambda)$. If we index such that $h_1$ is the first gridpoint, then the two BC's imply $h_{-1}=h_3$ and $h_0 = h_2$, and similarly at the other end of the domain. Can you tell me if the following, which is my approach, is correct?

I should say, since $h_0$ and $h_{-1}$ don't exist in our indexing, to satisfy the boundary conditions I concatenate $h$ by 4 components, 2 on each side. These are ghost points that I use to satisfy the BC. Code wrote in MATLAB.

However, a rough stability check implies a parameter $dt/dx^4$ which is very large for the specified $dx$ value. Does anyone have any ideas on how to handle this problem, perhaps a different way?

%% Parameters
hstar = 0.01;
n = 3;
m = 2;
theta = 50*pi/180;
K = 2*(1-cos(theta))/hstar;
G = 1;
eps = 10^(-3);
lambda = 2.66;

%% Numerics
tf = 200;% final time
dt = 0.01;% time step
t = 0:dt:tf;% time vector

xmin = 0;
xmax = 2.5;
dx = 0.01;% spatial step
x = xmin:dx:xmax;% x-vector

xc = xmin-2*dx:dx:xmax+2*dx;% concatenate h by 4 fictitious values (two at each
    % domain endpoint) to accommodate BC. Then h actually starts at i=3 and
    % ends at i = N-2.

%% IC
h0 = 0.1;% IC
hi = h0*(1+eps*cos(2*pi*xc/lambda));% IC
h = hi;

% preallocate
hnew = h;
h2x = zeros(size(x));
hxx = zeros(size(x));
hxxxx = zeros(size(x));
fh = zeros(size(x));
fhh = zeros(size(x));

N = length(h);

%% FTCS difference equations for eq. (4)
for j = 1:size(t,2)
    fh = (m*(hstar./h).^m-n*(hstar./h).^n)./h;% f'(h)
    fhh = ((hstar./h).^n*n*(n+1)-(hstar./h).^m*m*(m+1))./h.^2;% f''(h)
 
    % BC
    h(1) = h(5);
    h(2) = h(4);
 
    h(N) = h(N-4);
    h(N-1) = h(N-3);
     
    % FTCS
    for i = 3:N-2
        h2x(i) = (h(i+1)^2-h(i-1)^2)/(2*dx);% (h^2)'(x)
        hxx(i) = (h(i-1)-2*h(i)+h(i+1))/dx^2;% h''(x)
        hxxxx(i) = (h(i-2)-4*h(i-1)+6*h(i)-4*h(i+1)+h(i+2))/dx^4;% h''''(x)
     
        hnew(i) = h(i) - dt*h(i)^3*(hxxxx(i)+(K*fh(i)-G)*hxx(i)+K*h2x(i)*fhh(i));
    end% for i
 
    h = hnew;% update time
end% for j

h = h(3:N-2);% delete fictitious points used for BC

 

[member 428835]

Shoot, I can't edit the code above but here is an update:

%% Parameters
hstar = 0.01;
n = 3;
m = 2;
theta = 50*pi/180;
K = 2*(1-cos(theta))/hstar;
G = 1;
eps = 10^(-3);
lambda = 2.66;

%% Numerics

xmin = 0;
xmax = lambda;
dx = 0.01;% spatial step
x = xmin:dx:xmax;% x-vector

xc = xmin-2*dx:dx:xmax+2*dx;% concatenate h by 4 fictitious values (two at each
    % domain endpoint) to accommodate BC. Then h actually starts at i=3 and
    % ends at i = N-2.

%% IC
h0 = 0.1;% IC
hi = h0*(1+eps*cos(2*pi*xc/lambda));% IC
h = hi;

% preallocate
hnew = h;
hx = zeros(size(x));
h2x = zeros(size(x));
h3x = zeros(size(x));
hxx = zeros(size(x));
hxxx = zeros(size(x));
hxxxx = zeros(size(x));
fh = zeros(size(x));
fhh = zeros(size(x));
p1 = zeros(size(x));
p2 = zeros(size(x));

N = length(h);

dt = 0.01^4;% time step
tmax = 5;% max time
t = 0:dt:tmax;% time vector

%% FTCS difference equations for eq. (4)
for j = 1:size(t,2)
    fh = (m*(hstar./h).^m-n*(hstar./h).^n)./h;% f'(h)
    fhh = ((hstar./h).^n*n*(n+1)-(hstar./h).^m*m*(m+1))./h.^2;% f''(h)
  
    % BC
    h(1) = h(5);
    h(2) = h(4);
  
    h(N) = h(N-4);
    h(N-1) = h(N-3);
      
    % FTCS
    for i = 3:N-2
        hx(i) = (h(i+1)-h(i-1))/(2*dx);% h'(x)
        h2x(i) = (h(i+1)^2-h(i-1)^2)/(2*dx);% (h^2)'(x)
        h3x(i) = (h(i+1)^3-h(i-1)^3)/(2*dx);% (h^3)'(x)
        hxx(i) = (h(i-1)-2*h(i)+h(i+1))/dx^2;% h''(x)
        hxxx(i) = (-h(i-2)/2+h(i-1)-h(i+1)+h(i+2)/2)/dx^3;% h'''(x)
        hxxxx(i) = (h(i-2)-4*h(i-1)+6*h(i)-4*h(i+1)+h(i+2))/dx^4;% h''''(x)
      
        p1(i) = h3x(i)*(hxxx(i)+(K*fh(i)-G)*hx(i));% 1/2 product rule
        p2(i) = h(i)^3*(hxxxx(i)+K*fhh(i)*h2x(i)+(K*fh(i)-G)*hxx(i));% 2/2 product rule
      
        hnew(i) = h(i) - dt*(p1(i)+p2(i));
    end% for i
end% for j

h = h(3:N-2);% delete fictitious points used for BC

 

[hilbert2]

This seems to be the thin film equation with a disjoining pressure included... When I was studying that subject myself, I unfortunately didn't attempt to finite difference that equation myself but used a commercial finite element solver. There was some difficulty in forming a Crank-Nicolson scheme when there's that disjoining pressure that goes infinite at zero film thickness, if I remember correctly. That kind of equations are a quite specialized field, and it may be difficult to find anyone here who's familiar with them, but good luck.

 

[member 428835]

Thanks for your input hilbert2 (and you're spot on regarding the thin film)! Would you be able to tell me if you agree with the following discretization (where superscript j denotes time--only used on LHS for simplicity, as RHS is all at time j--and subscript i denotes space) for the following PDE, as this is really my question

$$
h_t+\partial_x\left[h^3(h_{xxx}+(Kf'(h)-G)h_x\right] = 0 \implies\\
h_t + \partial_x h^3 \left(h_{xxx}+(Kf'(h)-G)h_x\right) + h^3 \left(h_{xxxx}+Kf''(h)h_x^2+(Kf'(h)-G)h_{xx}\right) = 0 \implies\\
h^{j+1}_i = h_i - \Delta t \left[ \partial_x h^3 \left(h_{xxx}+(Kf'_i-G)h_x\right) + h^3 \left(h_{xxxx}+Kf_i''h_x^2+(Kf'(h)-G)h_{xx}\right) \right]:\\
h_x = \frac{h_{i+1}-h_{i-1}}{2\Delta x}\\
h_x^2 = \frac{h^2_{i+1}-h^2_{i-1}}{2\Delta x}\\
h_x^3 = \frac{h^3_{i+1}-h^3_{i-1}}{2\Delta x}\\
h_{xx} = \frac{h_{i+1}-2h_i+h_{i-1}}{\Delta x^2}\\
h_{xxx} = \frac{-h_{i-2}/2+h_{i-1}-h_{i+1}+h_{i+2}/2}{\Delta x^3}\\
h_{xxxx} = \frac{h_{i-2}-4h_{i-1}+6h_i-4h_{i+1}+h_{i+2}}{\Delta x^4}$$
Boundary conditions are $h_x=h_{xxx}=0$ at both ends of the $x$ domain (however long you want to make it) and some initial height, call it $h_0(1+\epsilon \cos(2 \pi x/\lambda)$. If we index such that $h_1$ is the first gridpoint, then the two BC's imply $h_{-1}=h_3$ and $h_0 = h_2$, and similarly at the other end of the domain. Can you tell me if the above is correct?