如何在MATLAB中实现1D Richards方程的Neumann-Dirichlet上边界切换

我在MATLAB中采用全隐式有限差分法（Picard迭代）求解一维混合形式Richards方程，z轴向下为正（z=0为土壤表面）。底部边界（z=80cm）设为恒定负压水头-51cm；顶部边界（z=0cm）为随时间变化的降雨通量q=-P(t)，且表面最大允许压头h=20cm。以下是我编写的处理该顶部边界条件的MATLAB代码：

% 1D Integer-order Richards' Equation (Picard Iteration Scheme)
clc;
clear; close all;

L = 80;                % Column length [cm]
T = 816;               % Total time [h]

Ks = 3.2;              % cm/h (Saturated hydraulic conductivity)
alpha_v = 0.076;       % 1/m  (van Genuchten parameter; if using cm throughout, strictly it should be converted to 1/cm)
nv = 1.86;
mv = 1 - 1/nv;
theta_s = 0.320;
theta_r = 0.032;

h_bot = -51;           % cm (Constant pressure head at the bottom boundary)

% Initial water content -> Initial pressure head
theta_init = 0.08;
Se = (theta_init-theta_r)/(theta_s-theta_r);
Psi_init = -( (Se^(-1/mv)-1)^(1/nv) )/alpha_v;  % cm (Negative)

dz = 0.1;              % cm
dt = 0.05;             % h

z = 0:dz:L;
t = 0:dt:T;
M = length(z);
N = length(t);

% Apply upper boundary pulse infiltration (Flux, positive downwards, unit: cm/h)
qTop = zeros(1,N);
q_pulse = 15/0.5;                  % 30 cm/h
events = [0 0.5; 240 240.5; 480 480.5];  % hours
for k = 1:size(events,1)
    qTop(t>=events(k,1) & t<=events(k,2)) = q_pulse;
end

% van Genuchten-Mualem Model
Kfun = @(h) (Ks .* ((1 - (alpha_v * abs(h)).^(nv-1) .* (1 + (alpha_v * abs(h)).^nv).^(-mv)).^2) ./ ...
               (1 + (alpha_v * abs(h)).^nv).^(mv/2));
Cfun = @(h) (theta_s - theta_r) .* alpha_v .* nv .* mv .* ((alpha_v .*abs(h)).^(nv-1)) .* ((1 + ( alpha_v *abs(h)).^nv).^(-mv-1));
Qfun = @(h) (theta_r + (theta_s - theta_r) .* (1 + (alpha_v *abs(h)).^nv).^(-mv));

% Result arrays
h = zeros(M, N);
theta = zeros(M, N);
q = zeros(M, N);

h(:,1) = Psi_init;
theta(:,1) = Qfun(Psi_init);

% Picard parameters
maxIter = 200;
tol = 1e-6;                  
iterate = zeros(N,1);

% Workspace for iterations
KK = zeros(M,1);
CC = zeros(M,1);

% Main time-stepping loop
for n = 1:N-1

    % Previous time step (converged solution)
    h_old = h(:,n);

    % Picard initial guess: use previous time step
    h_it = h_old;

    for it = 1:maxIter

        % Update nonlinear coefficients using current iteration solution (Picard: fix coefficients, solve linear equation for new h)
        for j = 1:M
            if h_it(j) > 0
                KK(j) = Ks;
                CC(j) = 0;
            else
                KK(j) = Kfun(h_it(j));
                CC(j) = Cfun(h_it(j));
            end
        end

        % Assemble tridiagonal linear system B*h_new = d
        B = zeros(M,M);
        d = zeros(M,1);

        % q = -K * (dh/dz + 1)
        h_max = 20;   % Maximum allowable pressure head (ponding depth)
        
        K12 = 0.5*(KK(1)+KK(2));
        qT = qTop(n+1);
        
          
         if qT > 0  % Water is being applied
            if h_it(1) < h_max
                B(1,1) = -1; B(1,2) = 1;
                d(1)   = dz * (qT/max(K12,1e-12) + 1); % [Correction] z is positive downwards, gravity term is -1
            else
                B(1,1) = 1; B(1,2) = 0; 
                d(1)   = h_max; % Ponding limit reached, force constant head
            end
        else       % Water application stopped
            if h_it(1) > 0
                B(1,1) = 1; B(1,2) = 0; 
                d(1)   = 0; % Surface is ponded, let it infiltrate under potential energy, set head to 0
            else
                B(1,1) = -1; B(1,2) = 1;
                d(1)   = dz * (0/max(K12,1e-12) + 1); % No ponding, zero flux boundary (qT=0)
            end
        end


        for j = 2:M-1
            Kjp = 0.5*(KK(j)+KK(j+1));   % K_{j+1/2}
            Kjm = 0.5*(KK(j)+KK(j-1));   % K_{j-1/2}

            B(j,j-1) = -(dt/(dz^2))*Kjm;
            B(j,j)   =  CC(j) + (dt/(dz^2))*(Kjm + Kjp);
            B(j,j+1) = -(dt/(dz^2))*Kjp;

            d(j) = CC(j)*h_old(j) - (dt/dz)*(Kjp - Kjm);
        end

        % Bottom boundary: Constant head
         B(M, M-1) = 0;
        B(M,M) = 1;
        d(M)   = h_bot;
        h_new = ThomasAlgorithm(B, d);

        % Convergence criterion
        err = max(abs(h_new - h_it));
        h_it = h_new;

        if err < tol
            iterate(n+1) = it;
            break;
        end
    end

    % Save converged solution for the current time step
    h(:,n+1) = h_it;

    % Calculate flux q
    for j = 1:M
        h_j = h(j, n+1);
        if h_j > 0
            K_j = Ks;
        else
            K_j = Kfun(h_j);
        end

        if j == 1
            dh_dz = (h(j+1, n+1) - h_j) / (-dz);
        elseif j == M
            dh_dz = (h_j - h(j-1, n+1)) / (-dz);
        else
            dh_dz = (h(j+1, n+1) - h(j-1, n+1)) / (-2*dz);
        end

        q(j, n+1) = -K_j * (dh_dz + 1);
    end
    
  for j = 1:M
       if h(j,n+1)<0
        theta(j,n+1)=Qfun(h(j,n+1));
       else
        theta(j,n+1)=theta_s;   
       
       end
  end

    fprintf('Picard-Richards: Time step %d / %d, Iterations: %d\n', n, N-1, max(1,iterate(n+1)));
end

% Thomas算法函数
function x = ThomasAlgorithm(A, d)
    % A 是三对角矩阵
    % d 是右端项
    n = length(d);
    
    % 创建向量
    c_prime = zeros(n, 1);
    d_prime = zeros(n, 1);
    
    % 前向消去
    c_prime(1) = A(1,2) / A(1,1);
    d_prime(1) = d(1) / A(1,1);
    
    for i = 2:n-1
        temp = A(i,i) - A(i,i-1) * c_prime(i-1);
        c_prime(i) = A(i,i+1) / temp;
        d_prime(i) = (d(i) - A(i,i-1) * d_prime(i-1)) / temp;
    end
    
    d_prime(n) = (d(n) - A(n,n-1) * d_prime(n-1)) / (A(n,n) - A(n,n-1) * c_prime(n-1));
    
    % 回代
    x = zeros(n, 1);
    x(n) = d_prime(n);
    
    for i = n-1:-1:1
        x(i) = d_prime(i) - c_prime(i) * x(i+1);
    end
end

内容的提问来源于stack exchange，提问作者user32543063