降雨条件下坡面径流和入渗耦合的数值模拟MATLAB 复现代码

2293801076 发代码源码

部分代码

for k=1:1:t  %按时间层循环
    %上游迁移法计算K(i+-1/2)、D(i+-1/2)
    K_p = K; K_m(1) = K(1); K_m(2:end) = K(1:end-1);
    D_p = D; D_m(1) = D(1); D_m(2:end) = D(1:end-1);
    
    %B,I,J向量赋值
    %B(i)为D_i = \theta_i^j-\frac{\Delta t}{\Delta z}[K(\theta)_{i+\frac{1}{2}}^{j+1}-K(\theta)_{i-\frac{1}{2}}^{j+1}]
    B = x - dz * lambda * (K_p - K_m);
    B(1) = x(1); B(N) = x(N);
    
    %三对角矩阵顶角向量I赋值，I(i)为方程组系数B(i)
    I = 1 + lambda * (D_p - D_m);
    I(1) = 1; I(N) = 1;
    
    %三对角矩阵上一向量J赋值，J(i)为方程组系数C(i)
    J = - lambda * D_p(1:end-1);
    J(1) = 0;
    
    %三对角矩阵下一向量M赋值，M(i)为方程组系数A(i)
    M = lambda * D_m(2:end);
    M(N-1) = 0;    
    
    %组合I,J,M向量得到三对角矩阵A
    A = diag(I) + diag(J,1) + diag(M,-1);
    
    %解五对角矩阵差分方程组，并以列形式存储在组合矩阵X中
    X(:,k+1) = A\B';
    
    %这一次的解作为下一次循环的初始值
    x=X(:,k+1)';
    
    %更新非饱和土壤导水率K(thtea)和扩散率D(theta)
    K = K_t * (x / theta_s).^m_k;
    D = a * exp(b * x) * 10^(-4) /60;
end