#math

Project

前言

學測前被人推坑過曲線插值，現在就來寫寫看。
被推坑的是多項式插值
但是這篇不是多項式，而是樣條函數。

介紹

跟多項式插值不一樣的點有

用低次的樣條函數比較不會產生高次多項式函數的”亂流”，且效果不錯

預測全球人口的多項式擬合，次數到8時會產生反常結果(source: 10-1 線性迴歸：曲線擬合)
並非一個函數用到底，而是利用取樣點之間作分割，分成小區間，利用某些特性做連結，以達到光滑的目的

上面是每個子區間都是次數=1的樣子，當然可以使用次數更高的函數

高中都學過「插值法」，用來算不在log表上的數字或是求多項式，而樣條插值是利用分段的函數，構成一個大函數，預測內插值。

線性樣條插值，最簡單的插值法，高中常用(source: 三次樣條插值)

線性插值雖然簡單，但無法滿足工程上「光滑」的需求，於是三次樣條就誕生了，最高次數為三能確保不過度扭曲，又可以確保光滑。

推導

假設有一複雜函數$f(x)$(難以評估)，需要用簡單函數$s(x)$去逼近，我們現在有$n+1$個取樣點$x_0, x_1$…$x_n$，由此分割出$n$個子區間，
$x_0$, $x_1$中間的樣條函數為$s_1(x)$，$x_{i-1}$, $x_i$中間的樣條函數為$s_i(x)$，$i \in [1,n]$
那該如何求出$s_i(x) = q_i x^3 + r_i x^2 + s_i x + t_i$ ?
我們有以下條件:

$s_i(x_i) = s_{i+1}(x_i)$
$s_i’(x_i) = s_{i+1}’(x_i)$
$s_i’’(x_i) = s_{i+1}’’(x_i)$
因為m=3，可知所有$s_i’’(x)$為同一直線，可寫出直線

$\frac{y-s_i’’(x_i)}{x-x_i}=\frac{s_i’’(x_{i-1})-s_i’’(x_i)}{x_{i-1}-x_i}$

整理一下

$s_i’’(x)=y=(s_i’’(x_{i-1})-s_i’’(x_i))\frac{x-x_i}{x_{i-1}-x_i}+s_i’’(x_i)$

積分，這邊將$s_i’’(x_i)$寫成$m_i$($m_i$是對應到$x_i$的$i$，因為所有的$s_i’’(x)$都是同一條直線，如果是對應到$s_i$，$m_i$就不會有變化)

$s_i’’(x) = (m_{i-1}-m_i)\frac{x-x_i}{x_{i-1}-x_i}+m_i$
$s_i’’(x) = \frac{x-x_i}{x_{i-1}-x_i}m_{i-1} + \frac{x_{i-1}-x}{x_{i-1}-x_i}m_i$
$s_i’(x) = \frac{(x-x_i)^2}{2(x_{i-1}-x_i)}m_{i-1}-\frac{(x_{i-1}-x)^2}{2(x_{i-1}-x_i)}m_i+C_1$
$s_i(x) = \frac{(x-x_i)^3}{6(x_{i-1}-x_i)}m_{i-1}+\frac{(x_{i-1}-x)^3}{6(x_{i-1}-x_i)}m_i+C_1x+C_2$

現在要求$C_1$, $C_2$，將$s_i(x_i)$寫成$y_i$(同樣，是對應$x_i$)，過程略…

$C_1 = \frac{y_i-y_{i-1}}{x_i-x_{i-1}}-\frac{x_i-x_{i-1}}{6}(m_i - m_{i-1})$
$C_2 = y_i-\frac{y_i-y_{i-1}}{x_i-x_{i-1}}x_i-\frac{x_{i-1}-x_i}{6}(m_ix_{i-1}-m_{i-1}x_i)$

利用算好的$C_1$, $C_2$將$s_i(x)$寫完整

$s_i(x)=\frac{(x-x_i)^3}{6(x_{i-1}-x_i)}m_{i-1}+\frac{(x_{i-1}-x)^3}{6(x_{i-1}-x_i)}m_i+ [\frac{y_i-y_{i-1}}{x_i-x_{i-1}}-\frac{x_i-x_{i-1}}{6}(m_i - m_{i-1})]x$
$+y_i-\frac{y_i-y_{i-1}}{x_i-x_{i-1}}x_i-\frac{x_{i-1}-x_i}{6}(m_ix_{i-1}-m_{i-1}x_i)$

現在只有$m_i$需要求解，這時需要用$s_i’(x_i) = s_{i+1}’(x_i)$，注意這時是要先把$s_i(x)$跟$s_{i+1}(x)$寫出來再代$x_i$，不然容易搞混。
整理等式極為繁瑣，在此直接給出最後整理好的結果

$\frac{6}{x_{i-1}-x_{i+1}}(\frac{y_i-y_{i-1}}{x_i-x_{i-1}}-\frac{y_{i+1}-y_i}{x_{i+1}-x_i}) = \frac{x_{i-1}-x_i}{x_{i-1}-x_{i+1}}m_{i-1} + 2m_i + \frac{x_i-x_{i+1}}{x_{i-1}-x{i+1}}m_{i+1}$

代換一下，等下比較好寫

$k_i = a_im_{i-1} + 2m_i + c_im_{i+1}$, $i \in [1, n-1]$ —(1)

我們有$m_i$, $i \in [0,n]$共$n+1$個未知數要求解，但是上面只有$n-1$個方程式，剩下兩個需要從邊界條件中得到
邊界條件分三種:

$s_1’(x_0)=f’(x_0)$, $s_n’(x_n)=f’(x_n)$
$s_1’’(x_0)=f’’(x_0)$, $s_n’’(x_n)=f’’(x_n)$
$s_1’(x_0)=s_n(x_n)$, $s_1’’(x_0)=s_n’’(x_n)$

邊界條件

邊界條件有許多種，每種都會對樣條函式產生不盡相同的影響。

source: 三次樣條插值

第一種:
先將$s_1’(x_0)$, $s_n’(x_n)$寫出來

$y_0’=f’(x_0)=s_1’(x_0)=\frac{x_0-x_1}{2}m_0-\frac{y_1-y_0}{x_1-x_0}-\frac{1}{6}(x_1-x_0)(m_1-m_0)$
$y_n’=f’(x_n)=s_n’(x_n)=-\frac{x_{n-1}-x_n}{2}m_n+\frac{y_n-y_{n-1}}{x_n-x_{n-1}}-\frac{1}{6}(x_n-x_{n-1})(m_n-m_{n-1})$

$y_0’+\frac{y_1-y_0}{x_1-x_0}=\frac{x_0-x_1}{2}m_0-\frac{1}{6}(x_1-x_0)(m_1-m_0)$
$y_n’-\frac{y_n-y_{n-1}}{x_n-x_{n-1}}=-\frac{x_{n-1}-x_n}{2}m_n-\frac{1}{6}(x_n-x_{n-1})(m_n-m_{n-1})$

$k_0=\frac{6}{x_0-x_1}(y_0’+\frac{y_1-y_0}{x_1-x_0})=2m_0+m_1$
$k_n=\frac{6}{x_{n-1}-x_n}(\frac{y_n-y_{n-1}}{x_n-x_{n-1}}-y_n’)=m_{n-1}+2m_n$

利用那$n-1$個條件與算好的邊界條件寫出矩陣
$\begin{bmatrix} 2 & 1 \\\ a_1 & 2 & c_1 & \\\ & a_2 & 2 & c_2 & \\\ & & & … \\\ & & & & & 1 & 2\end{bmatrix}\vec{m}=\vec{k}$
解矩陣的步驟就等下再寫。

第二種:

$s_1’’(x_0)=f’’(x_0)=m_0$
$s_n’’(x_n)=f’’(x_n)=m_n$

由上面的等式知

$k_1-a_1m_0=+2m_1+c_1m_2$
$k_{n-1}-c_{n-1}m_n=a_{n-1}m_{n-2}+2m_{n-1}$

寫出矩陣
$\begin{bmatrix} 2 & c_1 \\\ a_2 & 2 & c_2 & \\\ & a_3 & 2 & c_3 & \\\ & & & … \\\ & & & & & a_{n-1} & 2\end{bmatrix}\vec{m}= \begin{bmatrix}k_1-a_1f’’(x_0) \\\ k2 \\\ … \\\ k_{n-1}-c_{n-1}f’’(x_n) \end{bmatrix}$
注意這裡的$\vec{m}$是$m_1$到$m_{n-1}$。

第三種:

$s_1’(x_0)=s_n’(x_n)$
$s_1’’(x_0)=s_n’’(x_n)$

由第二個等式可知

$m_0 = m_n$

現在我們只有$n$個未知數需要求解
由第一個等式，知

$2(x_0-x_1+x_{n-1}-x_n)m_0+(x_0-x_1)m_1+(x_{n-1}-x_n)m_{n-1}=6(\frac{y_1-y_0}{x_1-x_0}+\frac{y_n-y_{n-1}}{x_n-x_{n-1}})$

代換

$a_0m_0+b_0m_1+c_0m_{n-1}=k_0$

再由(1)得

$k_{n-1}=a_{n-1}m_{n-2}+2m_{n-1}+c_{n-1}m_n$

寫出矩陣
$\begin{bmatrix} a_0 & b_0 & & … & & & c_0 \\\ a_1 & 2 & c_1 & \\\ & a_2 & 2 & c_2 & \\\ & & & … \\\ c_{n-1} & & & … & & a_{n-1} & 2\end{bmatrix}\vec{m}=\vec{k}$
$\vec{m}$是從$0~n-1$，$\vec{k}$同理。

解方程

算出來的矩陣為$A$，寫出等式$A\vec{m}=\vec{k}$
將$A$拆解成$LU$矩陣，
$L = \begin{bmatrix} d_0 & & & & & \\\ l_1 & d_1 & & \\\ & l_2 & d_2 & & \\\ & & & … \\\ & & & & l_n & d_n\end{bmatrix}$
$U = \begin{bmatrix} 1 & u_0 & & & \\\ & 1 & u_1 & \\\ & & 1 & u_2 & \\\ & & & … & & \\\ & & & & & u_{n-1} \\\ & & & & & 1 \end{bmatrix}$

根據使用的邊界條件不同而有所不同，而第一,二種步驟是相同的
先推出$d,l,u$的遞推關係，解出$\vec{w}$ ($U\vec{m}$)
在利用$U\vec{m}=\vec{w}$解出$\vec{m}$

第一種:

由$A=LU$可以觀察到

$d_0=2,u_0=\frac{1}{2}$
$l_i=a_i$
$d_i=2-u_{i-1}l_i$
$u_i=\frac{c_i}{d_i}$

由$L\vec{w}=\vec{k}$知

$w_0=\frac{k_0}{d_0}$
$w_i=\frac{k_i-l_iw_{i-1}}{d_i}$

再由$U\vec{m}=\vec{w}$得

$m_n=w_n$
$m_i=w_i-u_im_{i+1}$

第二種:

由$A=LU$可以觀察到

$d_1=2$
$l_i=a_i$
$d_i=2-l_iu_{i-1}$
$u_i=\frac{c_i}{d_i}$

由$L\vec{w}=\vec{k}$知

$w_1=\frac{k_1}{d_1}$
$w_i=\frac{k_i-l_iw_{i-1}}{d_i}$

再由$U\vec{m}=\vec{w}$得

$m_n=w_n$
$m_i=w_i-u_im_{i+1}$

第三種:

這個無法用上面的方法拆，只能用高斯$O(n^3)$求逆矩陣

現在有了具體的步驟就可以寫出程式了。
解第一二種矩陣的方法有個具體的名字: 追趕法，不多介紹，畢竟不是主角。

code:

(目前只有第一,二種邊界條件)

目標函數:

double f(double x)
{
    return exp(x);
}
double df(double x)
{
    double dx = 10e-8;
    if (isnan(f(x + dx)) || isinf(f(x + dx)))
        return (f(x) - f(x - dx)) / dx;
    return (f(x + dx) - f(x)) / dx;
}
double ddf(double x)
{
    double dx = 10e-8;
    if (isnan(df(x + dx)) || isinf(df(x + dx)))
        return (df(x) - df(x - dx)) / dx;
    return (df(x + dx) - df(x)) / dx;
}

採樣點, 邊界類型設定

constexpr double datSeed = 1, datStep = 0.1;
constexpr int dataCnt = 50;
constexpr double st = datSeed, ed = datSeed + ((double)dataCnt - 1.0) * datStep, step = 0.001;
constexpr int Type = 2;

void setSample(vector<double>& v, vector<double>& u)
{
    v.clear(); u.clear();

    double x = datSeed;
    for (int i = 0; i < dataCnt; x += datStep, ++i)
        v.emplace_back(x), u.emplace_back(f(x));
}

矩陣A, 向量k的設定

void bondType1_SetMat(vector<vector<double>>& A, vector<double>& k)
{
    k[0] = 6 * (SampleY[1] - SampleY[0] - df(SampleX[0])) / (g[1] * g[1]);
    k[SampleN] = 6 * ((SampleY[SampleN - 1] - SampleY[SampleN]) / g[SampleN] - df(SampleX[SampleN])) / g[SampleN];

    A.resize(SampleN + 1);
    for (int i = 0; i < SampleN + 1; ++i)
        A[i].resize(SampleN + 1);
    for (int i = 1; i < SampleN; ++i)
    {
        A[i][i - 1] = g[i] / (g[i] + g[i + 1]);
        A[i][i] = 2;
        A[i][i + 1] = g[i + 1] / (g[i] + g[i + 1]);
    }
    A[0][0] = 2, A[0][1] = 1, A[SampleN][SampleN] = 2, A[SampleN][SampleN - 1] = 1;
}
void bondType2_SetMat(vector<vector<double>>& A, vector<double>& k)
{
    k[1] -= (SampleX[0] - SampleX[1]) * ddf(SampleX[0]) / (SampleX[0] - SampleX[2]);
    k[SampleN - 1] -= (SampleX[SampleN - 1] - SampleX[SampleN]) * ddf(SampleX[SampleN]) / (SampleX[SampleN - 2] - SampleX[SampleN]);

    A.resize(SampleN + 1);
    for (int i = 0; i < SampleN + 1; ++i)
        A[i].resize(SampleN + 1);
    for (int i = 2; i < SampleN - 1; ++i)
    {
        A[i][i - 1] = g[i] / (g[i] + g[i + 1]);
        A[i][i] = 2;
        A[i][i + 1] = g[i + 1] / (g[i] + g[i + 1]);
    }
    A[1][1] = 2, A[1][2] = (SampleX[1] - SampleX[2]) / (SampleX[0] - SampleX[2]);
    A[SampleN - 1][SampleN - 1] = 2, A[SampleN - 1][SampleN - 2] = (SampleX[SampleN - 2] - SampleX[SampleN - 1]) / (SampleX[SampleN - 2] - SampleX[SampleN]);
}

解矩陣

void bondType1(vector<vector<double>>& A, vector<double>& m, vector<double>& k)
{
    vector<double> l, d, u, w;
    l.resize(SampleN + 1); d.resize(SampleN + 1); u.resize(SampleN + 1), w.resize(SampleN + 1);

    d[0] = A[0][0], u[0] = A[0][1] / d[0], w[0] = k[0] / d[0];
    for (int i = 1; i <= SampleN; ++i)
    {
        l[i] = A[i][i - 1];
        d[i] = A[i][i] - l[i] * u[i - 1];
        if (i < SampleN)
            u[i] = A[i][i + 1] / d[i];
        w[i] = (k[i] - l[i] * w[i - 1]) / d[i];
    }

    m[SampleN] = w[SampleN];
    for (int i = SampleN - 1; i >= 0; --i)
    {
        m[i] = w[i] - m[i + 1] * u[i];
    }
}
void bondType2(vector<vector<double>>& A, vector<double>& m, vector<double>& k)
{
    vector<double> l, d, u, w;
    l.resize(SampleN + 1); d.resize(SampleN + 1); u.resize(SampleN + 1), w.resize(SampleN + 1);

    d[1] = 2, u[1] = A[1][2] / d[1], w[1] = k[1] / d[1];
    for (int i = 2; i < SampleN; ++i)
    {
        l[i] = A[i][i - 1];
        d[i] = 2 - l[i] * u[i - 1];
        if (i < SampleN - 1)
            u[i] = A[i][i + 1] / d[i];
        w[i] = (k[i] - l[i] * w[i - 1]) / d[i];
    }

    m[SampleN - 1] = w[SampleN - 1];
    for (int i = SampleN - 2; i >= 1; --i)
    {
        m[i] = w[i] - m[i + 1] * u[i];
    }
    m[0] = ddf(SampleX[0]);
    m[SampleN] = ddf(SampleX[SampleN]);
}

樣條函數

double s(double x)
{
    if ((x < SampleX[0]) || (x > *SampleX.rbegin()))
        return 0.0;

    int i = 1;
    while (x > SampleX[i]) ++i;
    double res = 0.0;
    res = (x - SampleX[i]) * (x - SampleX[i]) * (x - SampleX[i]) * m[i - 1] + (SampleX[i - 1] - x) * (SampleX[i - 1] - x) * (SampleX[i - 1] - x) * m[i];
    res /= (6 * g[i]);
    res += ((SampleY[i - 1] - SampleY[i]) / g[i] + g[i] * (m[i] - m[i - 1]) / 6) * x;
    res += SampleY[i] - (g[i] * g[i] * m[i]) / 6 - ((SampleY[i - 1] - SampleY[i]) / g[i] + (m[i] - m[i - 1]) * g[i] / 6) * SampleX[i];

    return res;
}

主函式

int main()
{
    setSample(SampleX, SampleY);

    SampleN = SampleX.size() - 1;

    g.resize(SampleN + 1);
    for (int i = 1; i <= SampleN; ++i)
        g[i] = SampleX[i - 1] - SampleX[i];

    k.resize(SampleN + 1);
    for (int i = 1; i < SampleN; ++i)
        k[i] = ((SampleY[i - 1] - SampleY[i]) / g[i] + (SampleY[i + 1] - SampleY[i]) / g[i + 1]) * 6 / (g[i] + g[i + 1]);

    m.resize(SampleN + 1);

    if (Type == 1)
    {
        bondType1_SetMat(A, k);
        bondType1(A, m, k);
    }
    else if (Type == 2)
    {
        bondType2_SetMat(A, k);
        bondType2(A, m, k);
    }

    fstream fs;
    fs.open("function.txt", ios::out);

    for (double x = st; x <= ed; x += step)
    {
        fs << f(x) << endl;
        fs << s(x) << endl;
    }

    fs.close();

    cout << st << " " << ed << endl;

    return 0;
}

函式可視化

import numpy as np
import matplotlib.pyplot as plt

f = open("function.txt", 'r')

x = np.arange(1, 10.9, 0.001)

a = []
b = []

for cnt in x:
    s = f.readline()
    a.append(float(s.strip()))
    s = f.readline()
    b.append(float(s.strip()))

y1 = np.array(a) # real blue
#print(y1)
y2 = np.array(b) # predict orange
#print(y2)

plt.figure()
plt.xlabel("x")
plt.ylabel("y")

plt.plot(x, y1)
plt.plot(x, y2)