Crank-Nicolson格式求解二维抛物线方程

Mi@MI

已于 2023-12-05 10:50:50 修改

阅读量2.1k

点赞数 18

分类专栏：计算数学|数值解文章标签：算法线性代数数学建模 matlab

于 2023-11-25 18:11:10 首次发布

本文链接：https://blog.csdn.net/m0_67196062/article/details/134617200

版权

计算数学|数值解专栏收录该内容

7 篇文章 1 订阅

订阅专栏

前言

在初学数值解时候，卡住我最大的问题就是二维抛物线方程的数值求解，当时也在网上搜索了相关代码，但觉得都写的不太清楚（可能是我自己理解不够），于是我便以C语言代码为参考，自己写出了CN格式求解二维抛物线方程的MATLAB代码，出于时间原因，欧拉方法以及紧致差分甚至更高阶的方法，如有时间后续会继续更新。

抛物线方程的形式

经典的抛物线方程是基于热传导问题而给出来的，其二维的基本形式可以由下列偏微分方程给出：

有过一定数值计算基础的同学想必都知道CN格式在一维情形下的离散格式，类比于一维格式，同样可以推导出二维格式，至于具体的推导过程可以参考《微分方程数值解》一书，只不过形式更加复杂，求解困难，相比于一维格式，其系数矩阵不再是三对角矩阵而是块三对角矩阵，因此直接计算的话，求解困难并且计算量巨大，因此才有了ADI格式求解，其具体构造过程也可以参考上书，对于理论没那么好的同学或者刚刚入门的同学，建议从格式推导入手，然后进行编程，至于其稳定性和收敛性等证明不要求掌握，因此本文着重于代码编写而不侧重于理论分析。

交替方向隐（ADI）格式

ADI方法在上述书籍中主要有三种，但其构造思想都是一致的，都是转换成三对角矩阵，这样一来就可以采用成熟而高效的追赶法进行求解，在这里我给出Peaceman-Rachford格式的形式以及MATLAB代码，至于其具体推导过程仍可以参考上述书籍，其余格式可以参考此代码进行编写，读者可自行撰写，难度不会很大。

数值实验

数值算例如下：

主函数代码如下：

% CN格式求解二维抛物线方程
%空间以及时间的右端点
function [err]=CN_pw_two(hx,hy,tau)
Lx=1;Ly=1;T=1;
% 空间以及时间步长
% hx=1/60;hy=1/40;tau=1/20;
x=(0:hx:Lx)';y=(0:hy:Ly)';t=(0:tau:T)';
m=length(x);n=length(y);l=length(t);
r1=tau/(hx*hx);r2=tau/(hy*hy);
% 计算精确解
ue=zeros(m,n,l);
for i=1:m
    for j=1:n
        for k=1:l
            ue(i,j,k)=exp((x(i)+y(j))/2-t(k));
        end
    end
end
% 计算数值解
uc=zeros(m,n,l);
% 定义初始条件和边界条件
for i=1:m
    for j=1:n
        uc(i,j,1)=exp((x(i)+y(j))/2);
    end
end
for j=1:n
    for k=1:l
        uc(1,j,k)=exp(y(j)/2-t(k));
        uc(m,j,k)=exp((1+y(j))/2-t(k));
    end
end
for i=1:m
    for k=1:l
        uc(i,1,k)=exp(x(i)/2-t(k));
        uc(i,n,k)=exp((1+x(i))/2-t(k));
    end
end
% 构造系数矩阵
A=zeros(m-2);B=zeros(n-2);
for i=1:m-2
    A(i,i)=1+r1;
end
for i=1:m-3
    A(i,i+1)=-r1/2;A(i+1,i)=-r1/2;
end
for i=1:n-2
    B(i,i)=1+r2;
end
for i=1:n-3
    B(i,i+1)=-r2/2;B(i+1,i)=-r2/2;
end
d1=zeros(m-2,1);v=zeros(m-2,n-2);d2=zeros(n-2,1);
% 开始时间层循环,k
for k=1:l-1
% 求解V，固定j
    for j=1:n-2
        for i=1:m-2
            d1(i,1)=r2/2*(uc(i+1,j,k)+uc(i+1,j+2,k))+(1-r2)*uc(i+1,j+1,k)+tau/2*(-3/2)*(exp((x(i+1)+y(j+1))/2-(t(k)+t(k+1))/2));
        end
            v(1,j+1)=(1-r2)/2*uc(1,j+1,k)+(1+r2)/2*uc(1,j+1,k+1)+r2/4*(uc(1,j,k)+uc(1,j+2,k)-uc(1,j,k+1)-uc(1,j+2,k+1));
            v(m,j+1)=(1-r2)/2*uc(m,j+1,k)+(1+r2)/2*uc(m,j+1,k+1)+r2/4*(uc(m,j,k)+uc(m,j+2,k)-uc(m,j,k+1)-uc(m,j+2,k+1));
            d1(1,1)=d1(1,1)+r1/2*v(1,j+1);d1(end,1)=d1(end,1)+r1/2*v(m,j+1);           
            v1=sanduijiao(A,d1);
        for i=1:m-2
            v(i+1,j+1)=v1(i);
        end
    end
% 求解k+1时间层，固定i
    for i=1:m-2
        for j=1:n-2
            d2(j,1)=r1/2*(v(i,j+1)+v(i+2,j+1))+(1-r1)*v(i+1,j+1)+tau/2*(-3/2)*(exp((x(i+1)+y(j+1))/2-(t(k)+t(k+1))/2));
        end
            d2(1,1)=d2(1,1)+r2/2*uc(i+1,1,k+1);
            d2(end,1)=d2(end,1)+r2/2*uc(i+1,n,k+1);     
        v2=sanduijiao(B,d2);
        for j=1:n-2
            uc(i+1,j+1,k+1)=v2(j);
        end
    end
end
% 计算误差
err=abs(ue-uc);

收敛阶计算代码如下：

% 计算时间收敛阶
dt=[1/5,1/10,1/20,1/40,1/80];
l=length(dt);
for i=1:l
    et(i)=max(max(max(CN_pw_two(1/200,1/400,dt(i)))));
end
for i=1:l-1
    rate_t_CN_pw_two(i)=log2(et(i)/et(i+1));
end
% 计算空间收敛阶
dx=[1/5,1/10,1/20,1/40,1/80];
m=length(dx);
for i=1:m
    eh(i)=max(max(max(CN_pw_two(dx(i),1/500,1/2000))));
end
for i=1:m-1
    rate_h_CN_pw_two(i)=log2(eh(i)/eh(i+1));
end
% 由结果可以得到时间与空间x收敛阶确是二阶的，至于空间y也可以采用上述代码其结果仍是二阶，与理论分析一致

总结

由于上述代码中取得步长比较小，所以运行时间也会比较长，但是如果步长选取不当会导致收敛阶不对，其中的原因与理论分析有关，这里就不讲述吗，刚兴趣的伙伴们可以自行查找相关书籍分析其原因，至于图像的问题，本人执着于收敛阶的正确性，因此不过多关注于图像的展示，感兴趣的也可以自行编程画出图像。最后，本人所写代码追求简单而易懂，不过多参杂技巧性，如有不足之处望谅解。如对所有内容满意或者对数值解感兴趣的友友们，可以点个赞以及关注，后续我会继续更新相关的数值算法求解各种方程，你们的支持将是我持续更新的动力哒~