# A Finite Difference

Tags:
1. Oct 30, 2017

### joshmccraney

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?
Code (Text):

%% 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;% spacial 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 accomodate 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

Last edited: Oct 31, 2017
2. Oct 31, 2017

### joshmccraney

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

Code (Text):

%% 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;% spacial 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 accomodate 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

3. Nov 1, 2017

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

4. Nov 1, 2017

### joshmccraney

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?