有限元简述(3)：正交基函数与傅立叶级数

程序实例(I)

逼近函数 I 中介绍了逼近函数的数学推导, 这一次用程序实现它.

例题回顾

给定一个二次(抛物型)函数 $f(x) = 10(x-1)^2-1, x \in \Omega=[1, 2]$ , 在线性函数空间 $V$ 中寻找最佳逼近向量 $u$

V={\rm span}\{1, x\}.

$\bf 解:$ 设 $\psi_0 = 1, \psi_1 = x$ , 则

u = c_0\psi_0(x) + c_1\psi_1(x) = c_0 + c_1x.

所以系数矩阵

\begin{align*} & A_{11} = (\psi_0, \psi_0) = \int_1^2 1\cdot1{\rm d}x = 1,\\\\ & A_{12} = (\psi_1, \psi_0) = \int_1^2 x \cdot 1 {\rm d}x = \frac{3}{2},\\\\ & A_{21} = A_{12} = \frac{3}{2},\\\\ & A_{22} = (\psi_1, \psi_1) = \int_1^2 x \cdot x {\rm d}x = \frac{7}{3}. \end{align*}

所有系数矩阵是

A = \begin{bmatrix} 1 & \frac{3}{2} \\\\ \frac{3}{2} & \frac{7}{3} \end{bmatrix}

右端项为:

\begin{align*} & b_{1} = (\psi_0, f) = \int_1^2 1\cdot (10(x-1)^2-1){\rm d}x = \frac{7}{3},\\\\ & b_{2} = (\psi_1, f) = \int_1^2 x \cdot (10(x-1)^2-1) {\rm d}x =\frac{13}{3}. \end{align*}

所以 $b = [\frac{7}{3}, \frac{13}{3}]^{\sf T}$ , 解线性方程组就可以得到

c = \begin{bmatrix} −38/3 \\\\ 10\end{bmatrix}.

因此

u = 10x-38/3.

程序实现

主函数

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
close all
clear
clc
format rat % 以分式显示结果

% 基本变量
xx = [1, 2];
syms x
f = (10 * (x-1)^2 -1);
F = inline(vectorize(f), 'x');  % 符号向量化

% 基函数
N = 3;
for i = 1 : N
    V(i) = x^(i-1);
end

% 组装系数矩阵
A = assemble_stiffness_matrix(V, xx);

% 载入右端向量
b = assemble_load_vector(f, V, xx);

% 线性方程组求解
c = A\b;
u =V*c;

% 可视化
X = xx(1) : 0.01 : xx(2);
U = inline(vectorize(u), 'x');
Plot(X, F, U)

% 结果展示
result(A, b, c, u)

程序计算结果与手算的是一样的.

当 $V = {\rm span}\\\{1, x, x^2 \\\}$ 时, 逼近解等于真解, 如下图所示

本程序其他函数可点击原文链接下载, GitHub.

选取更好的基函数

上一节的基函数空间为 $V ={\rm span} \{ x^j\} , j\in {\mathcal{I}}_s,\, {\mathcal{I}}_s=\{0, 1, \dots, N \}$ , 在上一节的例子中, 函数逼近的很好, 在用此基函数逼近多项式时, 理论上可以得到原多项式. ..但是.., 当 $N$ 过大时, 形成的系数矩阵 $A$ 是奇异的, 是病态的, 即线性方程组系统不可解.

选择..正交..(或者几乎正交)的基函数是数值计算中经常使用的, 其原因是可以使得 $A_{ij}=0, i\neq j$ , 从而矩阵几乎是对角化的.

傅立叶级数

$\color{gray}{\textit{Fourier series}}$

令

V={\rm span} \{ \sin(\pi x), \sin(2\pi x), \dots, \sin(N+1)\pi x \}.

那么基函数为

\psi_i(x) = \sin(i+1)\pi x, \quad i \in \cal{I}_s.

将基函数带入上文的主程序中, 得到 N=3 和 N=11 时的拟合图

以上结果似乎拟合得很好, ..但是..可以发现, 无论当 $N$ 如何增大, 始终得到 $u(0)=u(1)=1$ . 肯定是哪里出错了:

u(x) = \sum_{j\in \mathcal{I}_s} c_j \sin(j+1)\pi x.

上式显示: $u(0) = u(1) \equiv 0$ . 因此需要修正算法:

令 $u(0)=f(0), u(1)=f(1)$ , 以加入边界信息, 再加上 $u(x) = \sum_{j\in \mathcal{I}_s} c_j \psi_j (x)$ , 可设

\tilde{u}(x) = (1-x)f(0) + xf(1) + \sum_{j\in \mathcal{I}_s} c_j \psi_j (x).

设 $B(x) = (1-x)f(0) + xf(1)$ , 此时的线性方程组系统为

\sum_{j\in \mathcal{I}_s} (\psi_j, \psi_i)c_j = (f-B, \psi_i), \quad i \in \mathcal{I}_s.

针对该基函数修正后的主函数为

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
close all
clear
clc
format rat % 以分式显示结果
% format long

% 基本变量
xx = [0, 1];
syms x
f = (10 * (x-1/2)^5 -1);
F = inline(vectorize(f), 'x');
f_0 = F(xx(1));
f_N = F(xx(end));

n = 6;
for i = 1:n
    V(i) = sin(i*pi*x);
end
B = f_0*(xx(end)-x) + f_N*(x-xx(1));

% 组装系数矩阵
A = assemble_stiffness_matrix(V, xx);

% 载入右端向量
b = assemble_load_vector(f-B, V, xx);

% 线性方程组求解
c = A\b;
u = B + V*c;

% 可视化
X = xx(1) : 0.01 : xx(2);
U = inline(vectorize(u), 'x');
Plot(X, F, U)

% 结果展示
result(A, b, c, u)