有限元入门： 用变分法求解一维常微分方程(一)

问题描述

给定下述方程
$$\{\begin{array}{rl}& -U^{\prime \prime }(t)+U(t)=({\pi }^2+1)\sin \pi t,\\ & U(0)=U(1)=0,\end{array}$$

其中 $t\in [0,1].$

• 请构造一次有限元求解上述方程并给出相应误差估计;
• 请构造二次有限元求解上述方程并给出相应误差估计.

问题分析

取检验函数空间 $V=H_0^1(\Omega ),$ 则对 $\forall v\in V,$ 有
$$\begin{array}{r}{\int }_{\Omega }(-u^{\prime \prime }+u)v\mathrm{d}t=({\pi }^2+1){\int }_{\Omega }\sin \pi t\cdot v\mathrm{d}t,\end{array}$$
由 Green 公式以及边值条件, 我们有
$${\int }_{\Omega }(\nabla u\nabla v+uv)\mathrm{d}t=({\pi }^2+1){\int }_{\Omega }\sin \pi t\cdot v\mathrm{d}t.$$
记 $$a(u,v)={\int }_{\Omega }(\nabla u\nabla v+uv)\mathrm{d}t,(f,v)=({\pi }^2+1){\int }_{\Omega }\sin \pi t\cdot v\mathrm{d}t.\text{\hspace{1em}(1)}$$
则方程(1) 可转化为
$$\{\begin{array}{rl}& Findu\in V,s.t.\\ & a(u,v)=(f,v),\forall v\in V.\end{array}\text{\hspace{1em}(2)}$$
设 $V_h$ 是 $V$ 的有限维逼近, 于是问题 (2) 可转化为
$$\{\begin{array}{rl}& Find{u}_h\in V_h,s.t.\\ & a(u_h,v_h)=(f,v_h)\forall v_h\in V_h.\end{array}\text{\hspace{1em}(3)}$$
取 $V_h$ 的一组基 { Φ i } i = 1 N , \{\Phi_i\}_{i=1}^N, {Φi​}i=1N​, 记 $u_h={\sum }_{i=1}^N{u}_i{\Phi }_i,$ 则 (3) 可转化为 $N$ 个方程
$$\begin{array}{r}a(\sum _{i=1}^N{u}_i{\Phi }_i,{\Phi }_j)=(f,{\Phi }_j),\forall j=1,2,\cdots ,N,\\ \sum _{i=1}^N{u}_{i}a({\Phi }_i,{\Phi }_j)=(f,{\Phi }_j),\forall j=1,2,\cdots ,N.\end{array}$$
令 $K_{ij}=a({\Phi }_i,{\Phi }_j),f_j=(f,{\Phi }_j),$ 则有
$$KU=F,$$
其中 $K=(k_{ij}{)}_{N\times N},F=(f_1,\cdots ,f_N{)}^T,U=(u_1,\cdots ,u_N{)}^T.$
记误差精度 $e(t)=U(t)-U_h(t),$ 其中 $U(t)$ 是精确解, $U_h(t)$ 是有限元解. 我们在此种模的意义下来度量误差
$$\Vert e{\Vert }^2={\int }_0^1|e^{\prime }(t){|}^2\mathrm{d}t.$$
由 \textcolor{red}{Ce$\acute{a}$ 引理} 可知, 存在常数 $C>0$, s.t.
$$\Vert e\Vert =\Vert U(t)-U_h(t)\Vert \le C\underset{v_h\in V_h}{\inf }\Vert U(t)-v_h(t)\Vert ,$$

一次有限元

数值格式

将区间 [0,1] 等分为 N 个网格, 记 $t_i=i/N,$ 故有 $t_0=0,t_N=1.$ 令 $h_i=t_i-t_{i-1}$ 为网格尺寸. 因为等分网格, 故 $h_i=h=1/N,\forall i.$ 取基函数 ${\Phi }_i(i=1,2,\cdots ,N-1)$ 为分片线性函数, 即
$${\Phi }_i(t)=\{\begin{array}{rlr}& \frac{t-t_{i-1}}{h},& t\in (t_{i-1},t_i],\\ & \frac{t_{i+1}-t}{h},& t\in (t_i,t_{i+1}],\\ & 0,& otherwise.\end{array}$$
令 [K_{ij}=\int_0^1 \Phi_i’(t)\Phi_j’(t)+\Phi_i(t)\Phi_j(t),dt,] 经计算有
$$K_{ij}=\{\begin{array}{rlr}& 0,& |i-j|>1,\\ & \frac{2}{h}+\frac{2}{3}h,& i=j,\\ & -\frac{1}{h}+\frac{1}{6}h,& otherwise.\end{array}$$
令 $$f_j={\int }_0^1({\pi }^2+1)\sin \pi t{\Phi }_j(t)\mathrm{d}t,$$
记 $u=[u_1,\cdots ,u_{N-1}{]}^T$, $F=[f_1,\cdots ,f_{N-1}{]}^T,$ 求解线性方程组
$$Ku=F,$$
便可用追赶法求得方程的数值解 $$U_h(t)=\sum _{i=1}^{N-1}{u}_i{\Phi }_i(t).$$
解常微分方程, 容易求得问题的解为
$$U(t)=\sin \pi t.$$
利用追赶法解线性方程组, 结果为

误差估计

仅需考虑插值误差估计即可. 事实上, 记$\omega (t)=U-\Pi U(t),$ 我们有
$$\begin{array}{rl}\Vert e{\Vert }^2& \le \Vert U-\Pi U(t){\Vert }^2={\int }_0^1|{\omega }^{\prime }(t){|}^2\mathrm{d}t=\sum _{i=1}^N{\int }_{t_{i-1}}^{t_i}|{\omega }^{\prime }(t){|}^2\mathrm{d}t\\ & =\sum _{i=1}^N{\int }_{t_{i-1}}^{t_i}|{\omega }^{\prime }(t)-{\omega }^{\prime }(\xi ){|}^2\mathrm{d}t=\sum _{i=1}^N{\int }_{t_{i-1}}^{t_i}({\int }_{\xi }^t|{\omega }^{\prime \prime }(s)|\mathrm{d}s{)}^2\mathrm{d}t\\ & =\sum _{i=1}^N{\int }_{t_{i-1}}^{t_i}({\int }_{\xi }^t|U^{\prime \prime }(s)|\mathrm{d}s{)}^2\mathrm{d}t\le \sum _{i=1}^N{\int }_{t_{i-1}}^{t_i}({\int }_{t_{i-1}}^{t_i}|U^{\prime \prime }(s)|\mathrm{d}s{)}^2\mathrm{d}t\\ & =\sum _{i=1}^N{h}_i({\int }_{t_{i-1}}^{t_i}|U^{\prime \prime }(s)|\mathrm{d}s{)}^2\\ & \le \sum _{i=1}^N{h}_i\left(\right({\int }_{t_{i-1}}^{t_i}|U^{\prime \prime }(\eta ){|}^2\mathrm{d}\eta {)}^{1/2}({\int }_{t_{i-1}}^{t_i}{1}^2\mathrm{d}\eta {)}^{1/2}{)}^2\\ & \le \sum _{i=1}^N{h}_i^2{\int }_{t_{i-1}}^{t_i}{U}^{\prime \prime }(\eta {)}^2\mathrm{d}\eta \le C{h}^2,\end{array}$$
即 $$\Vert e\Vert \le Ch,$$
具有一阶精度.
数值实验结果如下图所示, 使用 polyfit 拟合图像斜率有$$k=-0.9991,$$
即一次元的数值误差阶为 0.9991.

MATLAB程序

%% main program
n = 10;
err = zeros(8, 1);

% Linear interpolation
index = zeros(8, 1);
for i=1:8
    index(i) = 2^(i-1) * n;
end
for i = 1:8
    N = index(i);
    h = 1 / N;
    x = 0:h:1;
    xx = 0:1/1000:1;
    exact_solu = sin(pi.*xx);
    K = zeros(N-1);
    K = K + diag(ones(N-1,1).*(2/h+2/3*h), 0) + diag((h/6-1/h).*ones(N-2,1), 1) + diag((h/6-1/h).*ones(N-2,1), -1);
    F = (1+pi^2)/(h*pi^2) * (2.*sin(pi.*x(2:N))-sin(pi.*x(1:N-1))-sin(pi.*x(3:N+1)))';
    K = sparse(K);
    coeff = K\F;
    solu = [0; coeff; 0];
    figure()
    plot(x', solu, 'o--',xx', exact_solu, 'g--');
    legend('numerical solution', 'exact solution')
    xlabel('x')
    ylabel('U');
    xlim([-0.05,1.05])
    ylim([-0.05,1.05]);
    title(['numerical solution and exact solution N=',num2str(index(i))]);
 
% calculate the error
    du = solu(2:N+1) - solu(1:N);
    y1 = sin(2.*pi.*x);
    dy1 = y1(2:N+1) - y1(1:N);
    y2 = sin(pi.*x);
    dy2 = y2(2:N+1) - y2(1:N);
    error = pi^2*h/2*ones(1,N) + du'.^2/h + pi/4*dy1 - 2/h*dy2.*du';
    err(i) = sqrt(sum(error));
end
figure()
plot(log2(err),'o--');
title('误差阶')
xlabel('N')
ylabel('$log_{2}Error$', 'Interpreter','LaTex')
xlim([0.8,8.2])

axis on;
set(gca,'xtick',1:1:8);
set(gca,'xticklabel',{5,10,20,40,80,160,320,640});

[order,~] = polyfit(1:1:(length(err)), log2(err)', 1);
disp(order)

下一篇我们将介绍这个问题的二次元解法，详见：
https://blog.csdn.net/waitingwinter/article/details/106244868。