Hi PF!(adsbygoogle = window.adsbygoogle || []).push({});

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

**Physics Forums | Science Articles, Homework Help, Discussion**

Dismiss Notice

Join Physics Forums Today!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# A Finite Difference

Tags:

Have something to add?

Draft saved
Draft deleted

**Physics Forums | Science Articles, Homework Help, Discussion**