GMRES方法求解二维不可压Stokes方程(附MATLAB代码)

最新推荐文章于 2023-04-19 16:06:49 发布

waitingwinter

最新推荐文章于 2023-04-19 16:06:49 发布

阅读量6.6k

点赞数 6

分类专栏：计算数学文章标签： MATLAB GMRES 预优GMRES方法

本文链接：https://blog.csdn.net/waitingwinter/article/details/103863358

版权

计算数学专栏收录该内容

10 篇文章 14 订阅

订阅专栏

一、问题描述

考虑如下不可压的Stokes方程
$\left\{\begin{array}{cc} -\lambda \Delta \textbf{u} +\nabla p=\textbf{f}& in \;\Omega,\\ \nabla \textbf{u} = 0& in \; \Omega . \end{array}\right.$

通常设 $\textbf{f}=0,$ 令 $$\textbf{u}=(u,v)^T,\; \Gamma=\partial \Omega.$\\$
定义流函数 $$\psi$,$ 满足

$u=\frac{\partial \psi}{\partial y}, \; v=-\frac{\partial \psi}{\partial x}.$$

显然有 $\nabla \textbf{u}=0.$
定义窝度 $\omega$: $\omega = -\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x},$$ 经计算, stokes 方程可化为

$\left\{\begin{array}{cc} -\Delta u -\frac{\partial \omega}{\partial y}=0& in \;\Omega,\\ -\Delta v +\frac{\partial \omega }{\partial x}=0& in\;\Omega,\\ \lambda \Delta \omega =0& in\;\Omega. \end{array}\right.$

x相应的边界条件为

$\left\{\begin{array}{cc} u=1,\;v=0,\;\omega=-\frac{\partial u}{\partial y}& y=1,\; 0\leq x\leq 1,\\ u=0,\;v=0,\;\omega=-\frac{\partial u}{\partial y}& y=0,\; 0\leq x\leq 1,\\ u=0,\;v=0,\;\omega=\frac{\partial v}{\partial x}& x=1\;or\;0,\; 0\leq y\leq 1.\\ \end{array}\right.$

之后利用有限差分来计算.

问题：

1. 写出离散格式对应的代数系统, 即 $Ax=b$ ;

2. 设计有效的迭代法求解 $Ax=b,$ 要求单步迭代的复杂度为 $O(N),$ 其中 $N$ 为 $x$ 的维数;

3. 设计有效的预条件, 并比较预处理前后迭代法的效率差别.

二、问题分析

在[0,1]上进行等距差分, $\Delta x=\Delta y =1/(N-1).$ 记 $x_i=i\Delta x=i/(N-1),\;y_j=j\Delta y = j/(N-1),$ 其中 $$i,j=0,1,\cdots,N-1.$\\$
在 $(x_i,y_j)$ 点处, 有:

$\Delta u\approx \delta^2u_{ij}=\frac{u_{i+1,j}+u_{i-1,j}-2u_{ij}}{(\Delta x)^2}+\frac{u_{i,j+1}+u_{i,j-1}-2u_{i,j}}{(\Delta y)^2}, u_{i,j}\approx u(x_i,y_j),\\ \Delta v\approx \delta^2v_{ij}=\frac{v_{i+1,j}+v_{i-1,j}-2v_{ij}}{(\Delta x)^2}+\frac{v_{i,j+1}+v_{i,j-1}-2v_{i,j}}{(\Delta y)^2}, v_{i,j}\approx v(x_i,y_j),\\ \Delta \omega\approx \delta^2\omega_{ij}=\frac{\omega_{i+1,j}+\omega_{i-1,j}-2\omega_{ij}}{(\Delta x)^2}+\frac{\omega_{i,j+1}+\omega_{i,j-1}-2\omega_{i,j}}{(\Delta y)^2},\omega_{i,j}\approx \omega(x_i,y_j),$

$\frac{\partial \omega}{\partial x}\approx \delta _x\; \omega_{ij}=\frac{\omega_{i+1,j}-\omega_{i-1,j}}{2\Delta x},\quad \frac{\partial \omega}{\partial y}\approx \delta _y\; \omega_{ij}=\frac{\omega_{i,j+1}-\omega_{i,j-1}}{2\Delta y}.$

若 $(x_i,y_j)$ 是内点, 即 $1\leq i\leq N-2,\;1\leq j\leq N-2,$ 离散方程为

$\left\{\begin{array}{ccc} -\delta ^2 u_{ij}-\delta_y\; \omega_{ij}&=&0,\\ -\delta ^2 v_{ij}+\delta _x \;\omega_{ij}&=&0,\\ \delta ^2\;\omega_{ij}&=&0. \end{array}\right.$

若 $(x_i,y_j)$ 是边界点, 以 $i=0,1\leq j\leq N-2$ 为例

$\left\{\begin{array}{c} u_{0j}=0,\quad v_{0j}=0,\\ \omega_{0j}=-\frac{3v_{0,j}-4v_{1,j}+v_{2,j}}{2\Delta x}, \end{array}\right.$

其余边界点类似, 即

$\left\{\begin{array}{c} u_{N-1,j}=0,\;v_{N-1,j}=0,\quad j=1,\cdots,N-2.\\ \omega_{N-1,j}=\frac{3v_{N-1,j}-4v_{N-2,j}+v_{N-3,j}}{2\Delta x},\quad j=1,\cdots,N-2. \end{array}\right.$

$\left\{\begin{array}{c} u_{i,0}=0,\; v_{i,0}=0,\quad i=1,\cdots,N-2;\\ \omega_{i,0}=\frac{3u_{i,0}-4u_{i,1}+u_{i,2}}{2\Delta y} ,\quad \quad i=1,\cdots,N-2. \end{array}\right.$

$\left\{\begin{array}{c} u_{i,N-1}=1, \;v_{i,N-1}=0,\quad i=1,\cdots,N-2,\\ \omega_{i,N-1}=-\frac{3u_{i,N-1}-4u_{i,N-2}+u_{i,N-3}}{2\Delta y},\quad i=1,\cdots,N-2. \end{array}\right.$

对于角点

$\left\{\begin{array}{cc} u_{i,j}=1, & if \quad j=N-1,\\ u_{i,j}=v_{i,j}=\omega_{i,j}=0,& otherwise. \end{array}\right.$

三、数值格式

由第二部分的问题分析可知, 此问题对应的数值格式为 $Ax=b,$ 其中

$x=(u_{0,0},u_{1,0},\cdots,u_{N-1,N-1},v_{0,0},v_{1,0},\cdots,v_{N-1,N-1},\omega_{0,0},\omega_{1,0},\cdots,\omega_{N-1,N-1})^T,$

故可知 $x\in \mathbb{R}^{3N^2}.$

$A=\left[\begin{array}{ccc} A_{11}&A_{12}&A_{13}\\ A_{21}&A_{22}&A_{23}\\ A_{31}&A_{32}&A_{33} \end{array}\right],$

其中 $A_{ij}\in \mathbb{R}^{N^2}.$ 具体来说，考虑 $A_{11}$ 的元素, 经过分析易知, 这等价于 $-\Delta u$ 在区间 $[0,1]\times[0,1]$ 上的离散, 即有

$A_{11}=\left\{\begin{array}{cc} -1,& j=i,\;i=2,\cdots,N-1\;or \; i=N^2-N+2,\cdots,N^2-1;\\ -1,& j=i,\;i=(k-1)N+1,kN\;k=1,\cdots,N;\\ 4,& i=j,\; i=sN+t,\;s=2,\cdots,N-2,\;t=2,\cdots,N-1;\\ -1,& j=i\pm 1,\;i=sN+t,\;s=1,\cdots,N-2,\;t=2,\cdots,N-1;\\ -1,& j=i\pm N,\;i=sN+t,\;s=1,\cdots,N-2,\;t=2,\cdots,N-1;\\ 0,& otherwise. \end{array}\right.$

经分析有

$A_{22}=A_{33}=A_{11},\quad A_{12}=A_{21}=0.$

类似的办法来处理 $A_{13}$ 和 $A_{23},$ 有 $A_{13}$ 为 $-\delta _y\; \omega_{ij}$ 在 $[0,1]\times [0,1]$ 上的离散, $A_{23}$ 为 $\delta _x \;\omega_{ij}$ 在 $[0,1]\times [0,1]$ 上的离散. 即

$A_{13}=\left\{\begin{array}{cc} -\frac{h}{2},& j=i+N,\;i=sN+t,\;s=1,\cdots,N-2,\;t=2,\cdots,N-1;\\ \frac{h}{2},& j=i-N,\;i=sN+t,\;s=1,\cdots,N-2,\;t=2,\cdots,N-1;\\ 0,& otherwise. \end{array}\right.$

$A_{23}=\left\{\begin{array}{cc} \frac{h}{2},& j=i+1,\;i=sN+t,\;s=1,\cdots,N-2,\;t=2,\cdots,N-1;\\ -\frac{h}{2},& j=i-1,\;i=sN+t,\;s=1,\cdots,N-2,\;t=2,\cdots,N-1;\\ 0,&otherwise. \end{array}\right.$

$A_{31}=\left\{\begin{array}{cc} \frac{3}{2h},&j=i,\;i=N(N-1)+k,\;k=2,\cdots,N-1;\\ -\frac{2}{h},&j=i-N,\;i=N(N-1)+k,\;k=2,\cdots,N-1;\\ \frac{1}{2h},&j=i-2N,\;i=N(N-1)+k,\;k=2,\cdots,N-1;\\ -\frac{3}{2h},&j=i,\;i=k,\;k=2,\cdots,N-1;\\ \frac{2}{h},&j=i+N,\;i=k,\;k=2,\cdots,N-1;\\ -\frac{1}{2h},&j=i+2N,\;i=k,\;k=2,\cdots,N-1.\\ \end{array}\right.$

类似有

$A_{32}=\left\{\begin{array}{cc} \frac{3}{2h}&j=i,\;i=(k-1)N+1,k=2,\cdots,N-1;\\ -\frac{2}{h}&j=i+1,\;i=(k-1)N+1,k=2,\cdots,N-1;\\ \frac{1}{2h}&j=i+2,\;i=(k-1)N+1,k=2,\cdots,N-1;\\ -\frac{3}{2h}&j=i,\;i=kN,k=2,\cdots,N-1;\\ \frac{2}{h}&j=i-1,\;i=kN,k=2,\cdots,N-1;\\ -\frac{1}{2h}&j=i-2,\;i=kN,k=2,\cdots,N-1.\\ \end{array}\right.$
故我们得到系数矩阵 $A$ , 经分析易得 $b\in \mathbb{R}^{3N^2},$
$b(i)=\left\{\begin{array}{cc} 1,&N(N-1)+1<i<N^2;\\ 0,& otherwise. \end{array}\right.$

至此,我们便得到了离散的数值格式.

由图可知, $A$ 是不对称的稀疏矩阵, 故我们用GMRES算法和预优GMRES算法来求解此问题.

四、数值实验

GMRES方法

我们取[0,1]上的分隔点为5,10,20,40,80, 对应的矩阵 $A\in \mathbb{R}^{N\times N}$ 维数为 $$N=75,300,1200,4800,19200$.$

GMRES的MATLAB代码为

function [ x,j,res,resvec,time] = mygmres( A, b, x0, tol, max_it )
%% GMRES METHOD
% j : the number of iterations;
% res : the residual value at termination;
% x : solution
% time : CPU time;

if nargin < 5
    max_it=1000;
end
if nargin < 4
    tol = 1.e-10;
end
if nargin < 3
    x0 = zeros(length(b),1);
end
tic;
r0 = A*x0 - b;
beta = norm(r0);
V(:,1) = r0/beta;
resvec = beta;
Q = cast(1,'like',b);

j = 0;
while ((abs(resvec(j + 1)) > tol) && (j < max_it))
    % Arnoldi %
    j = j + 1;
    w = A*V(:,j);
    for i = 1:j
        H(i,j) = w'*V(:,i);
        w = w - H(i,j)*V(:,i);
    end
    H((j + 1),j) = norm(w,2);
    V(:,(j + 1)) = w/H((j + 1),j);
    
    % Construct R and Givens rotation %
    H(1:j,j) = Q*H(1:j,j);
    rho = H(j,j);
    H(j,j) = sqrt(rho^2 + H((j +1),j)^2);
    c = rho/H(j,j);
    s = H((j + 1),j)/H(j,j);
    H((j + 1),j) = 0;
    
    % Apply Givens rotation to Q %
    Q((j + 1),:) = -s*Q(j,:);
    Q(j,:) = c*Q(j,:);
    Q((j + 1),(j + 1)) = c;
    Q(j,(j + 1)) = s;
    
    % Apply Givens rotation to resvec%
    resvec(j + 1,1) = -s*resvec(j,1);
    resvec(j,1) = c*resvec(j,1);
end

y = H((1:j),:)\resvec(1:j);

x = x0 - V(:,(1:j))*y;
res = abs(resvec(j + 1));
resvec = abs(resvec);
time=toc;
end

预处理的GMRES方法

同样，取预处理矩阵为Jacobi矩阵,即

M = diag(diag(A));

有

比较两种方法, 我们有

其中, 预优GMRES方法为

function [ x,j,res,resvec,time] = mygmres_ja( A, b, x0, tol, max_it )
%% GMRES METHOD (Pre-processing M = Jacobi)
% j : the number of iterations;
% res : the residual value at termination;
% x : solution
% time : CPU time;

if nargin < 5
    max_it=1000;
end
if nargin < 4
    tol = 1.e-10;
end
if nargin < 3
    x0 = zeros(length(b),1);
end
tic;
M = diag(diag(A)); % Jacobi pre-processing matrix; 
r0 =M\(A*x0 - b);
beta = norm(r0);
V(:,1) = r0/beta;
resvec = beta;
Q = cast(1,'like',b);

j = 0;
while ((abs(resvec(j + 1)) > tol) && (j < max_it))
    % Arnoldi %
    j = j + 1;
    w = M\A*V(:,j);
    for i = 1:j
        H(i,j) = w'*V(:,i);
        w = w - H(i,j)*V(:,i);
    end
    H((j + 1),j) = norm(w,2);
    V(:,(j + 1)) = w/H((j + 1),j);
    
    % Construct R and Givens rotation %
    H(1:j,j) = Q*H(1:j,j);
    rho = H(j,j);
    H(j,j) = sqrt(rho^2 + H((j +1),j)^2);
    c = rho/H(j,j);
    s = H((j + 1),j)/H(j,j);
    H((j + 1),j) = 0;
    
    % Apply Givens rotation to Q %
    Q((j + 1),:) = -s*Q(j,:);
    Q(j,:) = c*Q(j,:);
    Q((j + 1),(j + 1)) = c;
    Q(j,(j + 1)) = s;
    
    % Apply Givens rotation to resvec%
    resvec(j + 1,1) = -s*resvec(j,1);
    resvec(j,1) = c*resvec(j,1);
end

y = H((1:j),:)\resvec(1:j);

x = x0 - V(:,(1:j))*y;
res = abs(resvec(j + 1));
resvec = abs(resvec);
time=toc;
end

数值解图像

最后, 我们给出最终的数值解图像

waitingwinter

关注

6
点赞
踩
40

收藏

觉得还不错? 一键收藏
2
评论
GMRES方法求解二维不可压Stokes方程(附MATLAB代码)

目录一、问题描述二、问题分析三、数值格式四、数值实验GMRES方法预处理的GMRES方法数值解图像一、问题描述考虑如下不可压的Stokes方程通常设令定义流函数满足显然有定义窝度经计算, stokes 方程可化为x相应的边界条件为之后利用有限差分来计算.问题：1. 写出离散格式对应的代数系...
复制链接

扫一扫