Peaceman-Rachford隐式格式求解二维抛物型PDE的MATLAB实现

对于抛物型偏微分的数值求解，常见的有向前/向后欧拉差分，Crank-Nicolson格式，peaceman-rachford格式等，对于二维偏微分方程的情形，比较好的策略是用P-R方法，其中用追赶法求解三对角矩阵。

虽然这种方法是比较好的，但是其数学形式较为复杂，并且多用Python语言来编写，在使用for循环时很容易出现索引的错误，导致编写的代码有问题。这里笔者给出使用MATLAB语言自编的P-R代码，其中大量使用矩阵运算，虽然比起Python可能更为晦涩难懂，但不至于迷惑于python的各种繁琐的数值运算操作，且运行速度大大提高，在缩短步长以达到更高精度的数值计算中可以起到比Python代码更为重要的作用

以下是代码实现，其中主函数是ADI_PR_method函数，只需要我们输入离散后的x,y,t向量即可，不过必须要求这些向量中的元素是等差数列，第一个元素是离散的起点，最后一个元素为离散的终点。返回值是一个三维的矩阵，每一个元素表示这个网格点处的近似值。

function result = ADI_PR_method(x,y,t)
%% x y t are arguements which is from the inf to the sup,we can create them by linspace
%% If you want to use this agorithm,you should replace the functions of u_xy0 in this m.file
%% as well as the latter five functions in peaceman_rachford.m file with your target functions
    x(1) = [];x(end) = [];y(1) = [];y(end) = [];
    hx = x(2)-x(1);hy = y(2)-y(1);tau = t(2)-t(1);
    rx = tau/hx^2;ry = tau/hy^2;
    u_k = u_xy0(x,y);
    result(:,:,1) = u_k;
    for k = 0:length(t)-1
        vij = zeros(length(x)+2,length(y));
        for j = 1:length(y)  
            vi = Thomas_method(x,y,j,rx,ry,k,tau,u_k,vij,'vi_compute');
            vij(:,j) = vi';
        end
        u_k = zeros(size(u_k,1),size(u_k,2));
        for i = 1:length(x)
            u_k(i+1,:) = Thomas_method(x,y,i,rx,ry,k,tau,u_k,vij,'uk_compute');
        end
        result(:,:,end+1) = u_k;
    end
end

function u_ij0 = u_xy0(x,y)
    x = [x(1)+x(1)-x(2),x,x(length(x))+x(2)-x(1)];
    y = [y(1)+y(1)-y(2),y,y(length(y))+y(2)-y(1)];
    [X,Y] = meshgrid(x,y);
    u_ij0 = exp((X+Y)/2);%you can only change this sentence in this function
end

function [a,b,c,f] = peaceman_rachford(x,y,num,rx,ry,k,tau,u_k,v_matrix,label)

    if label == 'vi_compute'
        a = ones(1,length(x)+2);a = a*(1+rx);a(1) = 1;a(end) = 1;
        b = ones(1,length(x)+1);b = b*(-rx/2);
        c = b;c(1) = 0;b(end) = 0;
        f = u_k(2:end-1,num:1:num+2)*[ry/2;1-ry;ry/2];
        f = f';
        f = f_xyt(num,x,y,k+1/2,tau,label)+f;

        f_1 = [(1-ry)/2,(1+ry)/2,ry/4]*[u_0yt(num,k,tau,y);u_0yt(num,k+1,tau,y);u_0yt(num-1,k,tau,y)+...
            u_0yt(num+1,k,tau,y)-u_0yt(num-1,k+1,tau,y)-u_0yt(num+1,k+1,tau,y)];

        f_2 = [(1-ry)/2,(1+ry)/2,ry/4]*[u_1yt(num,k,tau,y);u_1yt(num,k+1,tau,y);u_1yt(num-1,k,tau,y)+...
            u_1yt(num+1,k,tau,y)-u_1yt(num-1,k+1,tau,y)-u_1yt(num+1,k+1,tau,y)];
        f = [f_1,f,f_2]';
    elseif label == 'uk_compute'
        a = ones(1,length(y)+2);a = a*(1+ry);a(1) = 1;a(end) = 1;
        b = ones(1,length(y)+1);b = b*(-ry/2);
        c = b;c(1) = 0;b(end) = 0;
        f = [rx/2,1-rx,rx/2]*v_matrix(num:1:num+2,:);
        Fdt = f_xyt(num,x,y,k+1/2,tau,label);
        f = f+Fdt;
        f = [u_x0t(num,k+1,tau,x),f,u_x1t(num,k+1,tau,x)]';
    else
        fprint('输入的label格式错误');
    end
            
end


function f_ijk = f_xyt(num,x,y,k,tau,label)
    if label == 'vi_compute'
        t = k*tau;
        f_ijk = -3/2*exp((x+y(num))/2-t);%change here and you must obey the rule that x is a vector while y is a certain value in its vector
    elseif label == 'uk_compute'
        t = k*tau;
        f_ijk = -3/2*exp((x(num)+y)/2-t);%change here and you must obey the rule that y is a vector while x is a certain value in its vector
    end
    f_ijk = f_ijk*tau/2;% WARNING: this sentence can not be changed!
end

function u_0jk = u_0yt(j,k,tau,y) 
    if j==0
        yj = 0;
    elseif j==length(y)+1
        yj = y(end)+y(end)-y(end-1);
    else
        yj = y(j);
    end
    tk = k*tau;
    u_0jk = exp(yj/2-tk);%in this function,you can only change here
end

function u_1jk = u_1yt(j,k,tau,y)
    if j==0
        yj = 0;
    elseif j==length(y)+1
        yj = y(end)+y(end)-y(end-1);
    else
        yj = y(j);
    end
    tk = k*tau;
    u_1jk = exp((1+yj)/2-tk);%in this function,you can only change here
end

function u_i0k = u_x0t(i,k,tau,x)
    xi = x(i);
    tk = k*tau;
    u_i0k = exp(xi/2-tk);%in this function,you can only change here
end

function u_i1k = u_x1t(i,k,tau,x)
        xi = x(i);
    tk = k*tau;
    u_i1k = exp((1+xi)/2-tk);%in this function,you can only change here
end

function vector = Thomas_method(x,y,num,rx,ry,k,tau,u_k,v_matrix,label)
    [a,b,c,f] = peaceman_rachford(x,y,num,rx,ry,k,tau,u_k,v_matrix,label);
    n = length(a);
    q = zeros(1,n);q(1) = a(1);p = zeros(1,n-1);
    for flg = 2:n
        p(flg-1) = b(flg-1)/q(flg-1);
        q(flg) = a(flg)-p(flg-1)*c(flg-1);
    end
    temp = zeros(1,n);vector = temp;temp(1) = f(1);
    for flg = 2:n
        temp(flg) = f(flg)-p(flg-1)*temp(flg-1);
    end
    vector(end) = temp(end)/q(end);
    for flg = n-1:-1:1
        vector(flg) = (temp(flg)-c(flg)*vector(flg+1))/q(flg);
    end
end

本代码是以这个二维抛物型偏微分方程的数值求解为例来编写的，如果想移植本文的MATLAB代码用于你的实际问题的话，还需要对代码里的初值函数u_xy0以及Dirichlet边值函数u_0yt , u_1yt , u_x0t , u_x1t这些做出相应的修改。这里的命名方式相信读者能够看得懂，使用时请自行修改。

本文是参考华冬英老师的《微分方程的数值解法和程序实践》中的数学推导，并参考了Syphomn写的《二维空间的抛物型偏微分方程基本解法——ADI与紧ADI方法》写的Python程序，最终写出的，如有侵权，请联系我删除。

这是笔者的第一篇文章，各位读者有不同的意见或者想法，也欢迎评论或者私信我