Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

A Finite Difference

  1. Oct 30, 2017 #1

    joshmccraney

    User Avatar
    Gold Member

    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. jcsd
  3. Oct 31, 2017 #2

    joshmccraney

    User Avatar
    Gold Member

    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
     
     
  4. Nov 1, 2017 #3

    hilbert2

    User Avatar
    Science Advisor
    Gold Member

    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.
     
  5. Nov 1, 2017 #4

    joshmccraney

    User Avatar
    Gold Member

    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?
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted



Similar Discussions: Finite Difference
Loading...