strange behaviour in matlab while using finite difference method to solve pde

Hello everyone.

I am trying to solve the following pde

using the following boundary conditions at z = 0

the dot on top of phi_s means time derivative. The domain of the function is t = [0,160] seconds, z = [0,140] micrometers

I have tried implementing the finite difference method, but the solution is demonstrating weird behaviours, such as the one depicted below.

the value 0.0973 at phi_s(511,5001) jumped sharply to 0.9654 at phi_s(512,5001), although the derivative in the column direction at phi_s(511, 5001) is the same as the ones above it, so the change from phi_s(511,5001) to phi_s(512,5001) should be the same as the change from phi_s(510,5001) to phi_s(511,5001). I suspect that the value in phi_s(512,5001) is shifted by one decimal point by the program.

here is my finite difference code

clear
% parameters
K = 1e9;
G = 0.1 * K;
T = 160;
h_0 = 90e-6;
t = linspace(0, T, 1601);
z = linspace(0, h_0, 10001);
dt = t(2) - t(1);
dz = z(2) - z(1);
Nz = length(z);
Nt = length(t);

%initialize phi_s and dphi_s_dt
phi_s = zeros(Nt, Nz);
dphi_s_dt = zeros(Nt, Nz);

% define piecewise xi and eta
phi_cri = 0.03;
phi_w = 0.001;
eta_r = 1e11;
eta_g = 1.14*eta_r;
xi_r = 1;
xi_g = 5000*xi_r;
xi = @(phi_s) xi_r + 1/2 * (xi_g - xi_r) * (1 + tanh((phi_cri - phi_s)/phi_w));%(phi_s >= 0 & phi_s < 0.02) * xi_r + (phi_s >= 0.02) * xi_g;
eta = @(phi_s) eta_r + 1/2 * (eta_g - eta_r) * (1 + tanh((phi_cri - phi_s)/phi_w));%(phi_s >= 0 & phi_s < 0.02) * eta_g + (phi_s >= 0.02) * eta_r; % Piecewise eta

j = 1;
% define coefficients: A*u_xx + B*u_txx = C*u_t
A = (K + 4*G/3 * 1/0.68^2 * (0.68 - 0.32*0.36));
B = @(phi_s) eta(phi_s)/2 * (20/9 * (0.32^2/0.68^2 + 1) - 16/9*0.64/0.68^2);
C = @(phi_s) 2 * xi(phi_s) * 1.32^2;
D = 2 * K * 0.08 / 0.68; % Coefficient of dphi_s_dt
% solve z = 0 boundary condition
for k = 1:10
    for i = 1: Nt
        phi_s(i, 1) = D/A * (1 - exp(-A/B(phi_s(i,j)) * t(i)));
        dphi_s_dt(i, 1) = (D - A * phi_s(i, 1))/B(phi_s(i,j));
    end
end

% let u mean phi_s
u_xx = zeros(Nt, Nz - 2);
u_txx = zeros(Nt, Nz - 2);

% generating whole phi_s
for i = 1: Nt - 1
    if i == 1
        for j = 1: Nz - 2
            u_xx(i, j) = 1/dz^2 * (phi_s(i, j + 2) - 2*phi_s(i, j + 1) + phi_s(i, j));
            
        end
    end
    
    for j = 1: Nz - 2
        u_txx(i, j) = (C(phi_s(i,j)) * dphi_s_dt(i, j) - A * u_xx(i, j))/B(phi_s(i,j));
        u_xx(i + 1, j) = u_xx(i, j) + dt * u_txx(i, j);
        if j == 1
            u_x = (phi_s(i, j + 1) - phi_s(i, j))/dz;
            u_tx = -A/B(phi_s(i,j)) * u_x;
            dphi_s_dt(i, j + 1) = dz * u_tx + dphi_s_dt(i, j);
            phi_s(i + 1, j + 1) = dt * dphi_s_dt(i, j + 1) + phi_s(i, j + 1);
        end
        phi_s(i + 1, j + 2) = dz^2 * u_xx(i + 1, j) + 2 * phi_s(i + 1, j + 1) - phi_s(i + 1, j);
        phi_s(i + 1, j + 2) = max([0, min([0.2, phi_s(i + 1, j + 2)])]);
        dphi_s_dt(i, j + 2) = (phi_s(i + 1, j + 2) - phi_s(i, j + 2))/dt;
    end
end

Also the solution doesn't seem to converge no matter how high i increase the number of mesh points, in both t and z.

Is the sudden increase a bug in matlab or is it a bug in the code? How can I solve the convergence issue?

Thank you so much!

[John D'Errico]

NO. MATLAB does not have a bug in it as you are claiming, where it randomly/arbitrarily shifts numbers in an array by one decimal point, so by a power of 10? Even if that were conceivable, it would not happen, since MATLAB does not use base 10 to encode numbers anyway. And, given the uncounted numbers of times MATLAB has been used for the last 40 plus years, someone might have noticed a bug like that.

But your code may have a bug in it. That seems possible, in fact, even likely given the complexity of your code, and the fact you seem a relative novice to MATLAB programming.

You probably need to learn to use the debugger, looking carefully at those memory locations, and the things that might generate a discontinuity.

For example, a glance at your code shows some things that may be suspect. I see piecewise functions. Did you implement them properly? I see uses of max and min in your code, both things that explicitly create discontinuities. I see numbers that vary by many orders of magnitude in the coefficients, and then I see those numbers used in conjunction with functions like tanh, which scares me.

[Chin Ching]

Thank you for your comment!

Yes I've only been using matlab for 3 months, so im a complete novice. My code is bad for this reason.

Sadly the tanh and the cofficients that vary by orders of magnitude are needed. Maybe my method of solving the pde is already wrong, but i have no idea what is the most efficient way of solving this pde. May someone give me some insights?

Thanks a lot!

[Torsten]

Do you have a link to a publication where this kind of PDE is solved ?

Are you sure that both boundary conditions are to be imposed at z = 0 and no boundary condition at z = 140e-6 ?

What are your initial conditions ?

[Chin Ching]

Hi Torsten,

The boundary condition at z = 140e-6 is

The initial condition at is .

Here is a link to one of such publications: https://pubs.acs.org/doi/10.1021/acsmacrolett.4c00041

My equation is similar to equations 3 to 5 in this publication. The unknown functions there are u and α. . The generated graph of is Figure 2(a).

Thank you very much!

[S@m]

Hi Chin,

I ran your code in MATLAB R2024b in Windows machine. I did not observe the sudden jump you had mentioned (screenshot attached). Could you confirm if the issue is recurring at your end?