两点边值问题的有限元方法以及边值条件的处理示例

引言

本文参考李荣华教授的《偏微分方程数值解法》一书

题目

对于非齐次第二边值问题
$$\{\begin{array}{l}-\frac{d}{dx}(p\frac{du}{dx})+qu=f,a<x<b\\ u(a)=\alpha ,u^{\prime }(b)=\beta \end{array}$$
已知精确解$u(x)=\sin 5x+x^3(2-x)+2$，并且
$$\{\begin{array}{l}p(x)=\cos x\\ q(x)=x\end{array}$$
求剖分区间数分别为16和32时的精确解并画图

补全方程

使用$SageMath$进行符号运算，求出$u^{\prime }(x)$和$f(x)$，代码如下

x = var('x')
u = sin(5*x)+x^3*(2-x)+2
p = cos(x)
q = x
f = -derivative(p * derivative(u, x), x) + q * u
print(f)
print(derivative(u, x))

结果为：

-((x - 2)*x^3 - sin(5*x) - 2)*x + (6*(x - 2)*x + 6*x^2 + 25*sin(5*x))*cos(x) - (3*(x - 2)*x^2 + x^3 - 5*cos(5*x))*sin(x)
-3*(x - 2)*x^2 - x^3 + 5*cos(5*x)

于是原方程为
$$\{\begin{array}{ll}-\frac{d}{dx}(p\frac{du}{dx})+qu=-((x-2)x^3-\sin 5x-2)x+& \\ (6(x-2)x+6x^2+25\sin 5x)\cos (x)-(3(x-2)x^2+x^3-5\cos 5x)\sin x,& \\ & a<x<b\\ u(a)=\sin 5a+a^3(2-a)+2& \\ u^{\prime }(b)=-3(x-2)x^2-x^3+5\cos 5x& \end{array}$$

刚度矩阵

对于网格剖分
$$a=x_0<x_1<\cdots <x_n=b$$
对于每个单元$I_i=[x_{i-1},x_i]$，按照线性插值公式，有数值解
$$u_h=\frac{x_i-x}{h_i}{u}_{i-1}+\frac{x-x_{i-1}}{h_i}{u}_i,x\in I_i,i=1,2,\cdots ,n$$
其中$h_i=x_i-x_{i-1}$

为使按段插值标准化，通常用仿射变换
$$\xi =\frac{x-x_{i-1}}{h_i}$$
把$I_i$变到$\xi$轴上的参考单元$[0,1]$，令
$$N_0(\xi )=1-\xi ,N_i(\xi )=\xi$$
则
$$u_h=N_0(\xi )u_{i-1}+N_1(\xi )u_i$$
代入泛函
$$J(u)=\frac{1}{2}{\int }_a^b(p{u}^{\prime 2}+q{u}^2-2fu)dx$$
令
$$\frac{\partial J(u_h)}{\partial u_j}=a_{j-1,j}{u}_{j-1}+a_{jj}{u}_j+a_{j+1,j}{u}_{j+1}-bj=0$$
其中，$u_j=u_h(x_j),j=1,2,\cdots ,n,u_0=0$
$$\{\begin{array}{ll}{a}_{j-1,j}& ={\int }_0^1[-p(x_{j-1}+h_j\xi )/h_j+h_{j}q(x_{j-1}+h_j\xi )\xi (1-\xi )]d\xi \\ a_{j+1,j}& ={\int }_0^1[-p(x_j+h_{j+1}\xi )/h_{j+1}+h_{j+1}q(x_j+h_{j+1}\xi )\xi (1-\xi )]d\xi \\ a_{jj}& ={\int }_0^1[p(x_{j-1}+h_j\xi )/h_j+h_{j}q(x_{j-1}+h_j\xi ){\xi }^2]d\xi +\\ & {\int }_0^1[p(x_j+h_{j+1}\xi )/h_{j+1}+h_{j+1}q(x_j+h_{j+1}\xi )(1-\xi {)}^2]d\xi \\ b_j& ={\int }_0^1[h_{j}f(x_{j-1}+h_j\xi )\xi ]d\xi +{\int }_0^1[h_{j+1}f(x_j+h_{j+1}\xi )(1-\xi )]d\xi \end{array}\cdots (\ast )$$
于是可以得到总刚度矩阵$K$，满足
$$Ku=b$$
发现$a_{jj}$和$b_j$由不同的两项相加，于是可以得到单元刚度矩阵$K^{(i)}$
$$K^{(i)}=\left(\begin{array}{cc}{a}_{i-1,i-1}^{(i)}& a_{i-1,i}^{(i)}\\ a_{i,i-1}^{(i)}& a_{ii}^{(i)}\end{array}\right)$$
其中
$$\{\begin{array}{l}{a}_{i-1,i-1}^{(i)}=a_{i-1,i-1}^{(i)}={\int }_0^1[-p(x_{i-1}+h_i\xi )/h_i+h_{i}q(x_{i-1}+h_i\xi )\xi (1-\xi )]d\xi \\ a_{i-1,i-1}^{(i)}={\int }_0^1[p(x_{i-1}+h_i\xi )/h_i+h_{i}q(x_{i-1}+h_i\xi )(1-\xi {)}^2]d\xi \\ a_{ii}^{(i)}={\int }_0^1[p(x_i+h_{i+1}\xi )/h_{i+1}+h_{i+1}q(x_i+h_{i+1}\xi ){\xi }^2]d\xi \end{array}\cdots (\ast \ast )$$
不难发现，总刚度矩阵可以由单元刚度矩阵$K^{(i)}$“组装”而成，即对于刚度矩阵$K$中的元素$a_{ij}$和单元刚度矩阵$K^{(i)}$中的元素$a_{ij}^{(i)}$，有如下关系
$$a_{ij}=\{\begin{array}{ll}{a}_{i,i-1}^{(i)},& j=i-1\\ a_{ii}^{(i)}+a_{ii}^{(i+1)},& j=i\\ a_{i,i+1}^{(i+1)},& j=i+1\\ 0,& |j-i|\ge 2\end{array}$$
不难发现，总刚度矩阵$K$为三对角矩阵

对应的，向量$b$也可以进行“拆分”，令
$$\{\begin{array}{l}{f}_{i-1}^{(i)}={\int }_0^1[h_{j}f(x_{j-1}+h_j\xi )(1-\xi )]d\xi \\ f_i^{(i)}={\int }_0^1[h_{j+1}f(x_j+h_{j+1}\xi )\xi ]d\xi \end{array}$$
则有
$$b_1=f_1^{(1)}+f_1^{(2)},b_1=f_2^{(2)}+f_2^{(3)},\cdots ,b_1=f_n^{(n)}\cdots (\ast \ast \ast )$$

构造基底

再单元$I_i,I_{i+1}$考察线性插值公式及$u_i$的系数，我们自然对每一节点$x_i$构造山形函数：
$$\{\begin{array}{ll}{\phi }_i(x)& =\{\begin{array}{ll}1+\frac{x-x_i}{h_i},& x_{i-1}\le x\le x_i\\ 1-\frac{x-x_i}{h_{i+1}},& x_i\le x\le x_{i+1}\\ 0,& otherwise\end{array}\\ & i=1,2,\cdots ,n-1\\ {\phi }_n(x)& =\{\begin{array}{ll}1+\frac{x-x_n}{h_n},& x_{n-1}\le x\le x_n\\ 0,& otherwise\end{array}\end{array}$$
于是$u_h(x)=\sum _{i=1}^n{u}_i{\phi }_i(x),u_i=u_h(x_i)$

边值条件非齐次

左边值条件非齐次

$$u(\alpha )=a$$

则增加一基函数
$${\phi }_0(x)=\{\begin{array}{ll}a-\frac{x-x_0}{h_1},& x_0\le x\le x_1\\ 0,& otherwise\end{array}$$
将$u_h$表为
$$u_h=\alpha {\phi }_0(x)+\sum _{i=1}^n{u}_i{\phi }_i(x)$$
有限元方程为
$$\sum _{i=1}^{n}a({\phi }_i,{\phi }_j)u_i=(f,{\phi }_j)-\alpha a({\phi }_0,{\phi }_j),j=1,2,\cdots ,n$$
实际上只修改了第一个方程的右端，因为当$j=2,3,\cdots ,n$时，$a({\phi }_0,{\phi }_j)=0$

右边值条件非齐次

则右端修改为$(f,{\phi }_j)+p(b)\beta {\phi }_j(b)$，有限元方程变为
$$\sum _{i=1}^n{u}_i{\int }_a^b[p{\phi }_i^{\prime }{\phi }_j^{\prime }+q{\phi }_i{\phi }_j]dx=(f,{\phi }_j)+p(b)\beta {\phi }_j(b),j=1,2,\cdots ,n$$
实际上只修改了第$n$个方程的右端项，因为当$j=1,2,\cdots ,n-1$时，${\phi }_j(b)=0$

非齐次边值条件有限元方程

综上，我们得到非齐次边值条件有限元方程
$$\left(\begin{array}{ccc}a({\phi }_1,{\phi }_1)& \cdots & a({\phi }_n,{\phi }_1)\\ a({\phi }_1,{\phi }_2)& \cdots & a({\phi }_n,{\phi }_2)\\ ⋮& \ddots & ⋮\\ a({\phi }_1,{\phi }_{n-1})& \cdots & a({\phi }_n,{\phi }_{n-1})\\ a({\phi }_1,{\phi }_n)& \cdots & a({\phi }_n,{\phi }_n)\end{array}\right)\left(\begin{array}{c}{u}_1\\ u_2\\ ⋮\\ u_{n-1}\\ u_n\end{array}\right)=\left(\begin{array}{c}(f,{\phi }_1)-\alpha a({\phi }_0,{\phi }_1)\\ (f,{\phi }_2)\\ ⋮\\ (f,{\phi }_{n-1})\\ (f,{\phi }_n)+p(b)\beta {\phi }_n(b)\end{array}\right)$$

求数值解

我们以区间$[-1,2]$为例做均匀剖分

在求解数值解$u=(u_1,u_2,\cdots ,u_n{)}^T$时，我们可以直接求总刚度矩阵，也可以先求单元刚度矩阵再进行“组装”

直接求总刚度矩阵

$Matlab$代码如下

% 直接求刚度矩阵
n = ;
a = -1;
b_ = 2;

A = zeros(n, n);
b = zeros(n, 1);

p = zeros(1, n+1);
I = zeros(2, n);
for i=1:1:n+1
    p(1,i) = a + (b_-a)/n.*(i-1);
    if i<n+1
        I(1,i) = i;
        I(2,i) = i+1;
    end
end

u0 = @(x) sin(5.*x)+x.^3.*(2-x)+2;
du = @(x) -3.*(x - 2).*x.^2 - x.^3 + 5.*cos(5.*x);
alpha = u0(a);
beta = du(b_);
px = @(x) cos(x);
qx = @(x) x;
f = @(x) -((x - 2).*x.^3 - sin(5*x) - 2).*x + (6*(x - 2).*x + 6*x.^2 + 25*sin(5*x)).*cos(x) - (3*(x - 2).*x.^2 + x.^3 - 5*cos(5*x)).*sin(x);
aii = @(x1, x2, h) integral(@(ksi) (px(x1+h.*ksi)/h + h.*qx(x1+h.*ksi).*ksi.^2 + px(x2+h.*ksi)/h + h.*qx(x2+h.*ksi).*(1-ksi).^2), 0, 1);
a23 = @(x, h) integral(@(ksi) (-px(x+h.*ksi)/h + h.*qx(x+h.*ksi).*ksi.*(1-ksi)),0, 1);
ann = @(x, h) integral(@(ksi) (px(x+h.*ksi)/h + h.*qx(x+h.*ksi).*ksi.^2), 0, 1);
bi = @(x1, x2, h) integral(@(ksi) (h.*f(x1+h.*ksi).*ksi + h.*f(x2+h.*ksi).*(1-ksi)), 0, 1);
bn = @(x, h) integral(@(ksi) (h.*f(x+h.*ksi).*ksi), 0, 1);

for i=1:1:n-1
    x1 = p(1, I(1,i));
    x2 = p(1, I(2,i));
    h = x2 - x1;
    A(I(1,i), I(1,i)) = aii(x1, x2, h);
    b(I(1, i), 1) = bi(x1, x2, h);
end
for i=1:1:n-1
    x1 = p(1, I(1,i+1));
    x2 = p(1, I(2,i+1));
    h = x2 - x1;
    A(I(1,i), I(2,i)) = a23(x1, h);
    A(I(2,i), I(1,i)) = a23(x1, h);
end
x1 = p(1, I(1,n));
x2 = p(1, I(2,n));
h = x2 - x1;
A(n, n) = ann(x1, h);
% b(n, 1) = bi(x1, x2, h);
b(n, 1) = bn(p(1, I(1,n)), h);
% 边值条件
x1 = p(1, I(1,1));
b(I(1,1), 1) = b(I(1,1), 1) - alpha.*a23(x1, h);
b(I(2,n-1), 1) = b(I(2,n-1), 1) + px(b_).*beta;
u = A\b;
c = cond(A);

dx = 1/1000;
x = a:dx:b_;
u1 = zeros(1, length(x));
for i=1:1:n-1
    for j=1:1:length(x)
        if p(1, I(1,i))<=x(1, j) && x(1, j)<=p(1, I(2,i)) % x_{i-1} <= x <= x_i
            h = p(1, I(2, i)) - p(1, I(1, i));
            u1(1, j) = u1(1, j) + ((x(1, j) - p(1, I(1, i))) / h) * u(I(1, i), 1);
        elseif p(1, I(1,i+1))<=x(1, j) && x(1, j)<=p(1, I(2,i+1)) % x_i <= x <= x_{i+1}
            h = p(I(2, i+1)) - p(I(1, i+1));
            u1(1, j) = u1(1, j) + (1 - (x(1, j) - p(I(1, i+1))) / h) * u(I(1, i), 1);
        end
    end
end
for j=1:1:length(x)
    if p(1, I(1,n))<=x(1, j) && x(1, j)<=p(1, I(2,n)) % x_{n-1} <= x <= x_n
        h = p(1, I(2, n)) - p(1, I(1, n));
        u1(1, j) = u1(1, j) + ((x(1, j) - p(1, I(1, n))) / h) * u(I(1, n), 1);
    end
end
for j=1:1:length(x)
    if p(1, I(1,1))<=x(1, j) && x(1, j)<=p(1, I(2,1)) % x_0 <= x <= x_1
        h = p(1, I(2, n)) - p(1, I(1, n));
        u1(1, j) = u1(1, j) + alpha * (1 - (x(1, j) - p(1, I(1, 1))) / h);
    end
end

u_ = sin(5.*x)+x.^3.*(2-x)+2;

% 绘制数值解和精确解的图像
figure(1);
plot(x, u1, 'b', 'LineWidth', 1.5);
hold on;
plot(x, u_, 'g', 'LineWidth', 1.5);
xlabel('x');
ylabel('u(x)');
legend('数值解', '精确解');
title('两点边值问题的有限元方法');
grid on;
hold off;
box off;

当n=16时结果如下

当n=32时结果如下

先求单元刚度矩阵

$Matlab$代码如下

% 先求单元刚度矩阵
n = ;
a = -1;
b_ = 2;

A = zeros(n, n);
b = zeros(n, 1);

p = zeros(1, n+1);
I = zeros(2, n);
for i=1:1:n+1
    p(1,i) = a + (b_-a)/n.*(i-1);
    if i<n+1
        I(1,i) = i;
        I(2,i) = i+1;
    end
end

u0 = @(x) sin(5*x)+x.^3.*(2-x)+2;
du = @(x) -3*(x - 2).*x.^2 - x.^3 + 5*cos(5*x);
alpha = u0(a);
beta = du(b_);
px = @(x) cos(x);
qx = @(x) x;
f = @(x) -((x - 2).*x.^3 - sin(5*x) - 2).*x + (6*(x - 2).*x + 6*x.^2 + 25*sin(5*x)).*cos(x) - (3*(x - 2).*x.^2 + x.^3 - 5*cos(5*x)).*sin(x);
a1 = @(x, h) integral(@(ksi) (px(x+h.*ksi)/h + h.*qx(x+h.*ksi).*(1-ksi).^2), 0, 1);
a23 = @(x, h) integral(@(ksi) (-px(x+h.*ksi)/h + h.*qx(x+h.*ksi).*ksi.*(1-ksi)),0, 1);
a4 = @(x, h) integral(@(ksi) (px(x+h.*ksi)/h + h.*qx(x+h.*ksi).*ksi.^2), 0, 1);
b1 = @(x, h) integral(@(ksi) (h.*f(x+h.*ksi).*(1-ksi)), 0, 1);
b2 = @(x, h) integral(@(ksi) (h.*f(x+h.*ksi).*ksi), 0, 1);

for i=1:1:n-1
    x1 = p(1, I(1,i+1));
    x2 = p(1, I(2,i+1));
    h = x2 - x1;
    a1_ = a1(x1, h);
    a2_ = a23(x1, h);
    a3_ = a23(x1, h);
    a4_ = a4(x1, h);
    A(I(1,i), I(1,i)) = A(I(1,i), I(1,i)) + a1_;
    A(I(1,i), I(2,i)) = A(I(1,i), I(2,i)) + a2_;
    A(I(2,i), I(1,i)) = A(I(2,i), I(1,i)) + a3_;
    A(I(2,i), I(2,i)) = A(I(2,i), I(2,i)) + a4_;
end
x1 = p(1, I(1,1));
x2 = p(1, I(2,1));
h = x2 - x1;
A(I(1,1), I(1,1)) = A(I(1,1), I(1,1)) + a4(x1, h);
for i=2:1:n
    x1 = p(1, I(1,i));
    x2 = p(1, I(2,i));
    h = x2 - x1;
    b1_ = b1(x1, h);
    b2_ = b2(x1, h);
    b(I(1,i-1), 1) = b(I(1,i-1), 1) + b1_;
    b(I(2,i-1), 1) = b(I(2,i-1), 1) + b2_;
end
x1 = p(1, I(1,1));
x2 = p(1, I(2,1));
h = x2 - x1;
b(I(1,1), 1) = b(I(1,1), 1) + b2(x1, h);
% 边值条件
b(I(1,1), 1) = b(I(1,1), 1) - alpha.*a23(x1, h);
b(I(2,n-1), 1) = b(I(2,n-1), 1) + px(b_).*beta;
u = A\b;
c = cond(A);

dx = 1/1000;
x = a:dx:b_;
u1 = zeros(1, length(x));
for i=1:1:n-1
    for j=1:1:length(x)
        if p(1, I(1,i))<=x(1, j) && x(1, j)<=p(1, I(2,i)) % x_{i-1} <= x <= x_i
            h = p(1, I(2, i)) - p(1, I(1, i));
            u1(1, j) = u1(1, j) + ((x(1, j) - p(1, I(1, i))) / h) * u(I(1, i), 1);
        elseif p(1, I(1,i+1))<=x(1, j) && x(1, j)<=p(1, I(2,i+1)) % x_i <= x <= x_{i+1}
            h = p(I(2, i+1)) - p(I(1, i+1));
            u1(1, j) = u1(1, j) + (1 - (x(1, j) - p(I(1, i+1))) / h) * u(I(1, i), 1);
        end
    end
end
for j=1:1:length(x)
    if p(1, I(1,n))<=x(1, j) && x(1, j)<=p(1, I(2,n)) % x_{n-1} <= x <= x_n
        h = p(1, I(2, n)) - p(1, I(1, n));
        u1(1, j) = u1(1, j) + ((x(1, j) - p(1, I(1, n))) / h) * u(I(1, n), 1);
    end
end
for j=1:1:length(x)
    if p(1, I(1,1))<=x(1, j) && x(1, j)<=p(1, I(2,1)) % x_0 <= x <= x_1
        h = p(1, I(2, n)) - p(1, I(1, n));
        u1(1, j) = u1(1, j) + alpha * (1 - (x(1, j) - p(1, I(1, 1))) / h);
    end
end

u_ = sin(5.*x)+x.^3.*(2-x)+2;

% 绘制数值解和精确解的图像
figure(1);
plot(x, u1, 'r', 'LineWidth', 1.5);
hold on;
plot(x, u_, 'b', 'LineWidth', 1.5);
xlabel('x');
ylabel('u(x)');
legend('数值解', '精确解');
title('两点边值问题的有限元方法');
grid on;
hold off;
box off;

当n=16时结果如下

当n=32时结果如下