工程计算7——函数逼近与曲线拟合

连续函数空间C[a,b]

• 定义
  
  • 在区间[a,b]上连续函数的全体构成的集合
  • 空间是无限维的
  • 满足线性运算的封闭性
• 范数
  
  • $||f(x)||_\infty = max_{x\in[a,b]}|f(x)|$
  • $||f(x)||_1 = \int_a^b |f(x)|dx$
  • $||f(x)||_2 = [\int_a^b f^2(x)dx]^{1\over 2}$
• 内积
  
  • 离散：
    
    • $(f,g)=\sum_{i=0}^mw_if(x_i)g(x_i)$
    • $||f||_2 = \sqrt{(f,f)} = \sqrt{\sum_{i=0}^mw_i[f(x_i)]^2}$
  • 连续
    
    • $(f,g)=\int_a^b\rho(x_i)f(x_i)g(x_i)dx$
    • $||f||_2 = \sqrt{(f,f)} = \sqrt{\int_a^b\rho(x_i)[f(x_i)]^2dx}$

函数逼近

基本概念

• 定义
  
  • 在选定的一类函数中寻找某个函数g，使它是已知函数ƒ在一定意义下的近似表示，使得用g近似表示 ƒ而产生的误差最小。
• 分类
  
  • $\infty$范数：最佳一致逼近
  • 1范数、2范数
    
    • 曲线之间的面积最小
    • 2范数：最佳平方逼近和最小二乘逼近

连续函数的最佳平方逼近

• 定义
  
  • 
• 求解$S^*(x)$
  
  • 
  • 
  • 存在唯一性
• 最佳平方逼近多项式
  
  • 可能出现病态方程
  • 适用于低次函数逼近
• 用正交基求最佳平方逼近
  
  • 由于正交基的性质，G矩阵只有对角元不为0
  • $a_k^* = (\phi_k,f) / (\phi_k,\phi_k)$
  • 
  • 例：用勒让德多项式来求函数$f = \sqrt{1+x^2}$ $x\in [-1,1]$的逼近函数g

double g(int k,double v)
{
    if(k == 0) return 1;
    else if(k == 1) return v;
    else return (3 * v * v - 1) / 2.0;
}
double f(double x)
{
    return sqrt(x * x + 1);
}
double gf(int k)
{
    double s = 0;
    int nn = 10000;
    double begin_ = -1;
    for(int i = 0;i < nn; i++)
    {
        s += g(k,begin_) * f(begin_) * (2.0 / (nn * 1.0));
        begin_ += (2.0 / (nn * 1.0));
    }
    return s;
}
double gg(int k)
{
    double s = 0;
    int nn = 1000;
    double begin_ = -1;
    for(int i = 0;i < nn; i++)
    {
        s += g(k,begin_) * g(k,begin_)  * (2.0 / (nn * 1.0));
        begin_ += (2.0 / (nn * 1.0));
    }
    return s;
}
void Compute_AB()
{
    for(int i = 0; i < n; i++)
    {
        A[i] = gg(i);
    }
    for(int i = 0; i < n; i++)
    {
        B[i] = gf(i);
    }
}

曲线拟合的最小二乘法（离散）

• 定义
  
  • 已知函数值表$(x_i,y_i)$，在函数空间$\phi$求解$S^*(x)$，使得$\sum_{i=0}^mw_i[S^*(x_i) - y_i]^2 = min_{s(x)\in\phi}\sum_{i=0}^mw_i[S(x_i) - y_i]^2$
• 求解
  
  • 与连续函数得到一样的法方程
  • 只是$(\phi_k,f)$的计算方式不同。参考内积的离散和连续形式。
• 三项递推公式
  
  • 
• 用三项递推公式来拟合一下5个点

x	0	0.25	0.5	0.75	1
y	0.1	0.35	0.81	1.09	1.96

double f(double k,double v)
{
    if(k == 0) return 1;
    else if(k == 1) return v - 0.5;
    else return (v - 0.5) * (v - 0.5) - 0.125;
}
void Compute_B()
{
    double s;
    for(int i = 0; i < n; i++)
    {
        s = 0;
        for(int j = 0; j < m; j++)
        {
            s += f(i , x[j]) *  y[j];
        }
        B[i] = s;
    }
}
void Compute_A()
{
    double s;
    for(int i = 0; i < n; i++)
    {
        s = 0;
        for(int j = 0;j < m;j++)
        {
            s += f(i,x[j]) * f(i,x[j]);
        }
        A[i] = s;
    }
}