- #1
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?
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:
%% 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
Last edited by a moderator: