I have a PDE like this

$$ \frac{\partial h}{\partial t} = \bigg(\frac{\dot{L}}{L}\bigg)x\frac{\partial h}{\partial x} - \alpha\bigg[h^3\frac{\partial^3 h}{\partial x^3}\bigg] $$ With boundary and initial conditions given as

$$ h(0,t) = h(1,t) = 1; h_{xxx}(0,t) = 0; h_{xxx}(1,t) = \frac{1}{\alpha}\frac{\dot{L}}{L}; h(x,0) = 1$$

In this case, I have a uniform initial condition and when I use the ode15s solver to solve this PDE, the final result I get is always 1. So, the plot looks like the picture below but the original results are supposed to be curvy as time increases. Not sure what I am doing wrong and I hope someone can help. I have attached the code below.


clear all

pts = 2^10-1; tmax = 1.76;
ti = linspace(0,tmax,pts);
xi = linspace(0,1,pts);
h_init = ones(pts,1);

[t,y] = ode15s(@(t,y) example_model(t,y,xi'),ti,h_init);

tspan = [1 100 200 300 400 500];
for i = tspan
    hold on;

function dhdt = example_model(t,y,xi)
pts = length(y); h = y;
h([1 pts]) = 1; h0 = 100;
dx = xi(2)-xi(1);
% L and Lt
tau = 0.765; U0 = 0.37;
lambda = 11.6; L_cl = 0.4;
tt = t./tau; st = sqrt(tt);
term1 = -0.5.*(tt).^2;
term2 = 0.5.*sqrt(pi).*erf(st);
term3 = -st.*exp(-tt);
L = L_cl + U0.*tau.*(term1 + lambda.*(term2+term3));
Lt = -U0.*tau.*(t/tau^2 - lambda.*((exp(-tt).*st)./tau - exp(-tt)./(2.*tau.*st) + (3991211251234741.*exp(-tt))./(4503599627370496.*tau.*pi^(1/2).*st)));
if isnan(Lt)
    Lt = 0;
% constant
C = 7.87e-7; alpha = C*h0^3/(3*L^4);
% derivatives
hx = FD_derivative(h,dx);
hxx = FD_derivative(hx,dx);
hxxx = FD_derivative(hxx,dx);
q = (h.^3).*hxxx;
qx = FD_derivative(q,dx);
dhdx = (Lt/L).*xi.*hx;
dhdt = dhdx - alpha.*qx;

function dydx = FD_derivative(y,dx)
N = length(y); dydx = zeros(N,1);
for i = 1:N
    if i == 1
        dydx(i) = -y(i) + y(i+1);
    elseif i == N
        dydx(i) = y(i) - y(i-1);
        dydx(i) = -0.5*y(i-1) + 0.5*y(i+1);
dydx = dydx./dx;

1 Answer 1


Your PDE clearly states that if $h$ is uniform, it does not evolve (all spatial derivatives are zero). The only way to trigger an evolution is via the boundary conditions $h_{xxx}(0,t)$ and $h_{xxx}(1,t)$. It however seems that you never inject them in your duscretised PDE. Here your hxxx term at the boundaries does not involve the values of these imposed boundary conditions, since you only compute the gradients based in the inner points. My take would be to introduce ghost points (fake additional points, one on each side). Find the finite difference formula that uses them to compute $h_{XXX}$ on the sides, and use that to find what values of $h$ should be (dynamically) attributed to them so as to produce the desired values of $h_{XXX}$ at the sides.

  • 1
    $\begingroup$ Okay, this worked. The introduction of the hxxx boundary condition gives better solution $\endgroup$
    – dazemood
    Aug 2, 2022 at 7:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.