二阶椭圆型第一边值问题的数值解法（五点差分格式和有限体积法）附matlab代码

二阶椭圆型第一边值问题的数值解法（五点差分格式和有限体积法）附matlab代码

这里我们介绍五点差分格式和有限体积法求椭圆型第一边值问题，
题目：
分别采用矩形网格的五点差分格式和有限体积法求椭圆型第一边值问题，
$-\frac{{\partial }^{2}u}{\partial x^2}-\frac{{\partial }^{2}u}{\partial y^2}+a\frac{\partial u}{\partial x}=0,a>0,$
$u\left(0,y\right)=sin\left(\pi y\right),u\left(1,y\right)=0,0<y<1,$
$u\left(x,0\right)=u\left(x,1\right)=0,0<x<1.$

取参数 a = 40，得到不同网格最大误差和收敛阶。

问题分析

(i) 矩形网格的五点差分法:
Step1：网格剖分。取定沿 x,y 轴方向的步长分别为：$h_x,h_y$,便于计算我们取$h_x=h_y$,令 $h=\sqrt{h_x^2+h_y^2}$,作两族平行于 x,y 轴的直线$x=x_i,y=y_i$，其中$x_i=x_0+i{h}_x,i=0,1,2,\cdots ,N,y_j=y_0+j{h}_y,j=0,1,2,\cdots ,N$,两族直线的交点$(i{h}_x,j{h}_y)$称为网点或者节点,记为 $(x_i,y_j)$ .
Step2：微分算子离散.
在正则节点 $(x_i,y_j)$上沿 x,y 方向分别使用中心差商公式对 $u_{xx}$和 $u_{yy}$ 作近似，得到方程的五点差分格式为：
$-\left(\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h_x^2}+\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h_y^2}\right)+a\frac{u_{i+1,j}-u_{i-1,j}}{2h_x}=0$
$\Rightarrow$
$\left(-\frac{1}{h_x^2}-\frac{a}{2h_x}\right)u_{i-1,j}+\left(\frac{2}{h_x^2}+\frac{2}{h_y^2}\right)u_{i,j}+\left(-\frac{1}{h_x^2}+\frac{a}{2h_x}\right)u_{i+1,j}+\left(-\frac{1}{h_y^2}\right)u_{i,j-1}+\left(-\frac{1}{h_y^2}\right)u_{i,j+1}=0,i,j=1,2,\cdots ,N-1$
Let
A = [ 2 h x 2 + 2 h y 2 − 1 h x 2 + a 2 h x 0 ⋯ 0 0 − 1 h x 2 − a 2 h x 2 h x 2 + 2 h y 2 − 1 h x 2 + a 2 h x ⋯ 0 0 ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ 0 0 0 ⋯ − 1 h x 2 − a 2 h x 2 h x 2 + 2 h y 2 ] A=\left[\begin{matrix}\frac{2}{h_x^2}+\frac{2}{h_y^2}&-\frac{1}{h_x^2}+\frac{a}{2h_x}&0&\cdots&0&0\\-\frac{1}{h_x^2}-\frac{a}{2h_x}&\frac{2}{h_x^2}+\frac{2}{h_y^2}&-\frac{1}{h_x^2}+\frac{a}{2h_x}&\cdots&0&0\\\cdots&\cdots&\cdots&\cdots&\cdots&\cdots\\0&0&0&\cdots&-\frac{1}{h_x^2}-\frac{a}{2h_x}&\frac{2}{h_x^2}+\frac{2}{h_y^2}\\\end{matrix}\right] A= ​hx2​2​+hy2​2​−hx2​1​−2hx​a​⋯0​−hx2​1​+2hx​a​hx2​2​+hy2​2​⋯0​0−hx2​1​+2hx​a​⋯0​⋯⋯⋯⋯​00⋯−hx2​1​−2hx​a​​00⋯hx2​2​+hy2​2​​ ​
为一个$\left(N-1\right)\times \left(N-1\right)$矩阵,

B = [ − 1 h y 2 0 ⋯ 0 0 − 1 h y 2 ⋯ 0 ⋯ ⋯ ⋯ ⋯ 0 0 ⋯ − 1 h y 2 ] 是一个 ( N − 1 ) × ( N − 1 ) B=\left[\begin{matrix}-\frac{1}{h_y^2}&0&\cdots&0\\0&-\frac{1}{h_y^2}&\cdots&0\\\cdots&\cdots&\cdots&\cdots\\0&0&\cdots&-\frac{1}{h_y^2}\\\end{matrix}\right]是一个\left(N-1\right)\times\left(N-1\right) B= ​−hy2​1​0⋯0​0−hy2​1​⋯0​⋯⋯⋯⋯​00⋯−hy2​1​​ ​是一个(N−1)×(N−1)矩阵,
K = [ A B 0 ⋯ 0 0 0   B A B ⋯ 0 0 0 ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ 0 0 0 ⋯ 0 B A ] , K=\left[\begin{matrix}A&B&0&\cdots&0&0&0\\\ B&A&B&\cdots&0&0&0\\\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots\\0&0&0&\cdots&0&B&A\\\end{matrix}\right], K= ​A B⋯0​BA⋯0​0B⋯0​⋯⋯⋯⋯​00⋯0​00⋯B​00⋯A​ ​,是一个${\left(N-1\right)}^2\times {\left(N-1\right)}^2$矩阵,
$b_1={\left(\left(\frac{1}{h_x^2}+\frac{a}{2h_x}\right)u_{0,1}+{\frac{1}{h_y^2}u}_{1,0},{\frac{1}{h_y^2}u}_{2,0},\cdots ,{\frac{1}{h_y^2}u}_{N-2,0},{\frac{1}{h_y^2}u}_{N-1,0}+\left(\frac{1}{h_x^2}-\frac{a}{2h_x}\right)u_{N,1}\right)}^T$,
$b_i={\left(\left(\frac{1}{h_x^2}+\frac{a}{2h_x}\right)u_{0,i},0,\cdots ,0,\left(\frac{1}{h_x^2}-\frac{a}{2h_x}\right)u_{N,i}\right)}^T,i=2,3,\cdots ,N-2,b_{N-1}={\left(\left(\frac{1}{h_x^2}+\frac{a}{2h_x}\right)u_{0,N-1}+{\frac{1}{h_y^2}u}_{1,N},{\frac{1}{h_y^2}u}_{2,N},\cdots ,{\frac{1}{h_y^2}u}_{N-2,N},{\frac{1}{h_y^2}u}_{N-1,N}+\left(\frac{1}{h_x^2}-\frac{a}{2h_x}\right)u_{N,N-1}\right)}^T$,
$b={\left(b_1,b_2,\cdots ,b_{N-1}\right)}^T,U={\left(u_{1,1},u_{2,1},\cdots u_{N-1,1,}{u}_{1,2},u_{2,2},\cdots ,u_{N-1,2},\cdots u_{1,N-1},u_{2,N-1},\cdots u_{N-1,N-1}\right)}^T$.
因此，KU=b,解出U即可。
在本题中，因为
$u\left(0,y\right)=sin\left(\pi y\right),u\left(1,y\right)=0,0<y<1,u\left(x,0\right)=u\left(x,1\right)=0,0<x<1$.
所以
$b_1={\left(\left(\frac{1}{h_x^2}+\frac{a}{2h_x}\right)sin\left(\pi y_1\right),0,\cdots ,0,0\right)}^T,b_i={\left(\left(\frac{1}{h_x^2}+\frac{a}{2h_x}\right)sin\left(\pi y_i\right),0,\cdots ,0,0\right)}^T,i=2,3,\cdots ,N-2,b_{N-1}={\left(\left(\frac{1}{h_x^2}+\frac{a}{2h_x}\right)sin\left(\pi y_{N-1}\right),0,\cdots ,0,0\right)}^T$.

(ii) 有限体积法
我们需要作对偶剖分。记$x_{i-\frac{1}{2}}=\left(i-\frac{1}{2}\right)h_x,y_{j-\frac{1}{2}}=\left(j-\frac{1}{2}\right)h_y$
对于任一正则内点$(x_i,y_j)$ . 考虑对偶剖分的网点：$A(x_{i-\frac{1}{2}},y_{j-\frac{1}{2}}),B(x_{i+\frac{1}{2}},y_{j-\frac{1}{2}}),C(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}}),D(x_{i-\frac{1}{2}},y_{j+\frac{1}{2}})$,用$G_{ij}$表示以A,B,C,D为顶点的矩形域，$\partial G_{ij}$为其边界，利用格林公式，我们有
$-{\int }_{\partial G_{ij}}\frac{\partial u}{\partial n}+\int {\int }_{G_{ij}}a\frac{\partial u}{\partial x}dxdy=0$
$\Rightarrow$
$\frac{u_{i,j}{-u}_{i,j-1}}{h_y}{h}_x+\frac{u_{i,j}-u_{i+1,j}}{h_x}{h}_y+\frac{u_{i,j}{-u}_{i,j+1}}{h_y}{h}_x+\frac{u_{i,j}{-u}_{i-1,j}}{h_x}{h}_y+\frac{a}{2}{(u}_{i+1,j}-u_{i-1,j})h_y=0$
$\Rightarrow$
$(-{\frac{h_y}{h_x}-\frac{a}{2}{h}_y)u}_{i-1,j}+{\left(\frac{2h_x}{h_y}+\frac{2h_y}{h_x}\right)u}_{i,j}+(-{\frac{h_y}{h_x}+\frac{a}{2}{h}_y)u}_{i+1,j}-\frac{h_x}{h_y}{u}_{i,j-1}-\frac{h_x}{h_y}{u}_{i,j+1}=0$
Let
A = [ 2 h x h y + 2 h y h x − h y h x + a 2 h y 0 ⋯ 0 0 − h y h x − a 2 h y 2 h x h y + 2 h y h x − h y h x + a 2 h y ⋯ 0 0 ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ 0 0 0 ⋯ − h y h x − a 2 h y 2 h x h y + 2 h y h x ] A=\left[\begin{matrix}\frac{2h_x}{h_y}+\frac{2h_y}{h_x}&-\frac{h_y}{h_x}+\frac{a}{2}h_y&0&\cdots&0&0\\-\frac{h_y}{h_x}-\frac{a}{2}h_y&\frac{2h_x}{h_y}+\frac{2h_y}{h_x}&-\frac{h_y}{h_x}+\frac{a}{2}h_y&\cdots&0&0\\\cdots&\cdots&\cdots&\cdots&\cdots&\cdots\\0&0&0&\cdots&-\frac{h_y}{h_x}-\frac{a}{2}h_y&\frac{2h_x}{h_y}+\frac{2h_y}{h_x}\\\end{matrix}\right] A= ​hy​2hx​​+hx​2hy​​−hx​hy​​−2a​hy​⋯0​−hx​hy​​+2a​hy​hy​2hx​​+hx​2hy​​⋯0​0−hx​hy​​+2a​hy​⋯0​⋯⋯⋯⋯​00⋯−hx​hy​​−2a​hy​​00⋯hy​2hx​​+hx​2hy​​​ ​
为一个$\left(N-1\right)\times \left(N-1\right)$矩阵,
B = [ − h x h y 0 ⋯ 0 0 − h x h y ⋯ 0 ⋯ ⋯ ⋯ ⋯ 0 0 ⋯ − h x h y ] 为一个 ( N − 1 ) × ( N − 1 ) B=\left[\begin{matrix}-\frac{h_x}{h_y}&0&\cdots&0\\0&-\frac{h_x}{h_y}&\cdots&0\\\cdots&\cdots&\cdots&\cdots\\0&0&\cdots&-\frac{h_x}{h_y}\\\end{matrix}\right] 为一个\left(N-1\right)\times\left(N-1\right) B= ​−hy​hx​​0⋯0​0−hy​hx​​⋯0​⋯⋯⋯⋯​00⋯−hy​hx​​​ ​为一个(N−1)×(N−1)矩阵,
K = [ A B 0 ⋯ 0 0 0   B A B ⋯ 0 0 0 ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ 0 0 0 ⋯ 0 B A ] , K=\left[\begin{matrix}A&B&0&\cdots&0&0&0\\\ B&A&B&\cdots&0&0&0\\\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots\\0&0&0&\cdots&0&B&A\\\end{matrix}\right], K= ​A B⋯0​BA⋯0​0B⋯0​⋯⋯⋯⋯​00⋯0​00⋯B​00⋯A​ ​, 是一个${\left(N-1\right)}^2\times {\left(N-1\right)}^2$矩阵,
$b_1={\left((\frac{a}{2}{h}_y+\frac{h_y}{h_x}{)u}_{0,1}+{\frac{h_x}{h_y}u}_{1,0},{\frac{h_x}{h_y}u}_{2,0},\cdots ,{\frac{h_x}{h_y}u}_{N-2,0},{\frac{h_x}{h_y}u}_{N-1,0}+(\frac{h_y}{h_x}-\frac{a}{2}{h}_y)u_{N,1}\right)}^T,$
$b_i={\left((\frac{a}{2}{h}_y+\frac{h_y}{h_x}{)u}_{0,i},0,\cdots ,0,(\frac{h_y}{h_x}-\frac{a}{2}{h}_y)u_{N,i}\right)}^T,i=2,3,\cdots ,N-2,b_{N-1}={\left((\frac{a}{2}{h}_y+\frac{h_y}{h_x}{)u}_{0,N-1}+\frac{h_x}{h_y}{u}_{1,N},\cdots ,{\frac{h_x}{h_y}u}_{N-2,N},{\frac{h_x}{h_y}u}_{N-1,N}+(\frac{h_y}{h_x}-\frac{a}{2}{h}_y)u_{N,N-1}\right)}^T,$
$b={\left(b_1,b_2,\cdots ,b_{N-1}\right)}^T,U={\left(u_{1,1},u_{2,1},\cdots u_{N-1,1,}{u}_{1,2},u_{2,2},\cdots ,u_{N-1,2},\cdots u_{1,N-1},u_{2,N-1},\cdots u_{N-1,N-1}\right)}^T.$

因此，KU=b,解出U即可。

在本题中，因为
$u\left(0,y\right)=sin\left(\pi y\right),u\left(1,y\right)=0,0<y<1,u\left(x,0\right)=u\left(x,1\right)=0,0<x<1.$
所以
$$\begin{gathered} b_1={\left(\left(\frac{a}{2}{h}_y+\frac{h_y}{h_x}\right)sin\left(\pi y_1\right),0,\cdots ,0,0\right)}^T, \\ {b}_i={\left(\left(\frac{a}{2}{h}_y+\frac{h_y}{h_x}\right)sin\left(\pi y_i\right),0,\cdots ,0,0\right)}^T,i=2,3,\cdots ,N-2, \\ {b}_{N-1}={\left(\left(\frac{a}{2}{h}_y+\frac{h_y}{h_x}\right)sin\left(\pi y_{N-1}\right),0,\cdots ,0,0\right)}^T. \end{gathered}$$

为了得到两种格式的收敛阶，我们分别取N=8,16,32,64, 求出不同的N对应的最大误差
${ϵ}_{max}^N=max\left(\left|v_{i,j}-u_{i,j}\right|\right),i,j=1,2,\cdots ,N-1.其中，u_{i,j}为(x_i,y_j)处的数值解，v_{i,j}为(x_i,y_j)处的解析解，由于分割份数N的增量均为2倍关系，即以2为底的N次幂增加，由数值分析理论知识可知，收敛阶计算公式为：s_N=lo{g}_2\left(\frac{{ϵ}_{max}^N}{{ϵ}_{max}^{2N}}\right)$

全部代码如下(matlab)

function [E,U,xx,yy]=Fivepoints(N)%五点差分法
fun=inline('sin(pi*y)');
x0=0;
x1=1
y0=0;
y1=1;
h=1/N;
x=zeros(1,N);
x(1)=x0+h;
for i=1:(N-1)
    x(i+1)=x(i)+h;
end
y=zeros(1,N);
y(1)=y0+h;
for i=1:(N-1)
    y(i+1)=y(i)+h;
end
A=zeros(N-1,N-1);
A(1,1)=4/(h*h);
A(1,2)=-1/(h*h)+20/h;
A(N-1,N-1)=4/(h*h);
A(N-1,N-2)=-1/(h*h)-20/h;
for i=2:(N-2)
        A(i,i)=4/(h*h);
        A(i,i+1)=-1/(h*h)+20/h;
        A(i,i-1)=-1/(h*h)-20/h;
end
B=zeros(N-1,N-1);
for i=1:(N-1)
    B(i,i)=-1/(h*h);
end
K=kron(diag(ones(1,N-2),1)+diag(ones(1,N-2),-1),B)+kron(diag(ones(1,N-1)),A);
b=zeros(1,(N-1)*(N-1));
for i=1:(N-1)
    b((i-1)*(N-1)+1)=(1/(h*h)+20/h)*feval(fun,y(i));
end
U=K\b';
for i=1:(N-1)
    for j=1:(N-1)
        yy(1,(i-1)*(N-1)+j)=y(i);
    end
end
for i=1:(N-1)
    for j=1:(N-1)
        xx(1,(i-1)*(N-1)+j)=x(j);
    end
end
fun1=inline('sin(pi.*y).*(-exp(20+sqrt(400+pi*pi))/(exp(20-sqrt(400+pi*pi))-exp(20+sqrt(400+pi*pi)))*(exp(20-sqrt(400+pi*pi))).^x+exp(20-sqrt(400+pi*pi))/(exp(20-sqrt(400+pi*pi))-exp(20+sqrt(400+pi)))*(exp(20+sqrt(400+pi))).^x)','x','y');
for i=1:(N-1)*(N-1)
    z(i)=feval(fun1,xx(i),yy(i));
    e(i)=abs(U(i)-z(i));
end
E=max(e);

function [E,U,xx,yy]=V2(N)%有限体积法
fun=inline('sin(pi*y)');
x0=0;
x1=1
y0=0;
y1=1;
h=1/N;
x=zeros(1,N);
x(1)=x0+h;
for i=1:(N-1)
    x(i+1)=x(i)+h;
end
y=zeros(1,N);
y(1)=y0+h;
for i=1:(N-1)
    y(i+1)=y(i)+h;
end
A=zeros(N-1,N-1);
A(1,1)=4;
A(1,2)=-1+20*h;
A(N-1,N-1)=4;
A(N-1,N-2)=-1-20*h;
for i=2:(N-2)
        A(i,i)=4;
        A(i,i+1)=-1+20*h;
        A(i,i-1)=-1-20*h;
end
B=zeros(N-1,N-1);
for i=1:(N-1)
    B(i,i)=-1;
end
K=kron(diag(ones(1,N-2),1)+diag(ones(1,N-2),-1),B)+kron(diag(ones(1,N-1)),A);
b=zeros(1,(N-1)*(N-1));
for i=1:(N-1)
    b((i-1)*(N-1)+1)=(20*h+1)*feval(fun,y(i));
end
U=K\b';
for i=1:(N-1)
    for j=1:(N-1)
        yy(1,(i-1)*(N-1)+j)=y(i);
    end
end
for i=1:(N-1)
    for j=1:(N-1)
        xx(1,(i-1)*(N-1)+j)=x(j);
    end
end
fun1=inline('sin(pi.*y).*(-exp(20+sqrt(400+pi*pi))/(exp(20-sqrt(400+pi*pi))-exp(20+sqrt(400+pi*pi)))*(exp(20-sqrt(400+pi*pi))).^x+exp(20-sqrt(400+pi*pi))/(exp(20-sqrt(400+pi*pi))-exp(20+sqrt(400+pi)))*(exp(20+sqrt(400+pi))).^x)','x','y');
for i=1:(N-1)*(N-1)
    z(i)=feval(fun1,xx(i),yy(i));
    e(i)=abs(U(i)-z(i));
end
E=max(e);

%调试
clear all;clc;close  all;
[E,U1,xx1,yy1]=Fivepoints(20);%五点差分法
[E,U2,xx2,yy2]=V2(20);%有限体积法
[X1,Y1,Z1]=griddata(xx1,yy1,U1,linspace(0,1,20)',linspace(0,1,20),'v4');%插值
[X2,Y2,Z2]=griddata(xx2,yy2,U2,linspace(0,1,20)',linspace(0,1,20),'v4');%插值
subplot(2,2,1);
surf(X1,Y1,Z1)%三维曲面
title('五点差分法');
hold on;
subplot(2,2,2);
surf(X2,Y2,Z2)%三维曲面
title('有限体积法');
[x,y]=meshgrid(0:1/20:1);
z=sin(pi.*y).*(-exp(20+sqrt(400+pi*pi))/(exp(20-sqrt(400+pi*pi))-exp(20+sqrt(400+pi*pi)))*(exp(20-sqrt(400+pi*pi))).^x+exp(20-sqrt(400+pi*pi))/(exp(20-sqrt(400+pi*pi))-exp(20+sqrt(400+pi)))*(exp(20+sqrt(400+pi))).^x);
subplot(2,2,3);
mesh(x,y,z);%画出精确解
title('解析解');
N=[8,16,32,64];
for i=1:4
    e1(i,1)=Fivepoints(N(i));
    e2(i,1)=V2(N(i));
end
for i=1:3
    s1(i,1)=log2(e1(i,1)/e1(i+1,1));
    s2(i,1)=log2(e2(i,1)/e2(i+1,1));
end

问题结果

(N=20,h=0.05时，两种方法与解析解的图像如下所示）

从图中我们发现，数值解与解析解的图像高度吻合。为了得到两种方法的收敛阶，我们求出N=8,16,32,64对应的最大误差，并用收敛阶计算公式：$s_N=lo{g}_2\left(\frac{{ϵ}_{max}^N}{{ϵ}_{max}^{2N}}\right)$计算收敛阶。


从第一张表中我们可以看出，两种方法的最大误差高度相似，初步判定两种方法拥有一样的收敛阶。接着利用收敛阶计算公式，由计算结果可见，随着N值的增大，两种方法的收敛阶越来越接近于2，从而说明了五点差分法与有限体积法的收敛阶为2阶。文章来源地址https://www.toymoban.com/news/detail-421426.html

到了这里，关于二阶椭圆型第一边值问题的数值解法（五点差分格式和有限体积法）附matlab代码的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！