逼近函数 II

程序实例(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. $$

程序实现

主函数

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=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. $$

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

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)