前言
傅裡葉變換是貫穿一系列課程的應用，從起初的不知其是以然到學習電路理論時的輪廓逐漸清晰。
在電路理論的結課複習階段，我在圖書館找到國外教材翻閱，關注到行文中提到的吉布斯現象（Gibbs Phenomenon）【幸運發現了維基百科的鏡像網站，注意收藏】
這一小小的問題給我留下了深刻的印象，下決心對此鑽研。
過程
從錢學森開始
工程力學的課程資料中附有錢學森關于工程科學家的講稿，使我獲益匪淺，對我的個人發展有高度的指引意義，下面是他對工程科學家培養課程的見解：
隻能遺憾的說，目前的我，在很多方面尚未達到正常工程師的要求，但後面的話也稍稍使我寬慰：
顯 然 ， 對 于 一 個 有 發 展 前 途 的 工 程 科 學 家 來 說 ， 不 能 希 望 把 它 的 學 習 都 擠 在 四 年 大 學 裡 ． 事 實 上 ， 他 必 須 在 高 中 畢 業 後 進 入 一 個 好 的 工 學 院 ： 先 用 三 年 時 間 學 習 一 般 的 工 科 課 目 ： 然 後 他 必 須 用 大 約 三 年 時 間 學 習 科 學 和 數 學 ． 這 樣 一 來 ， 在 高 中 畢 業 後 至 少 要 用 六 年 時 間 來 培 訓 一 名 工 程 科 學 家 ； 而 現 在 的 實 際 情 況 是 ： 常 規 培 訓 隻 用 四 年 時 間
得益于綜合型考察的聯考制度，他提到的第一點我尚可符合，當然我也深知自己數理基礎在一定程度上的薄弱。當我不可免俗得尋求任何一個表現自己機會的同時，心中也有一個聲音一直提醒我——要注重終身發展。這段話也再次警醒我：職業規劃當長遠，數理基礎要夯實。
吉布斯相率與吉布斯現象
吉布斯相率
講稿中提到了冶金工藝的‘吉布斯相率’，再次喚起了我的回憶，在錢老的激勵下，我決定不放過這次研究機會，把二者搞清楚！
百度百科中關于吉布斯相率的解釋：
隻看公式并沒有清晰的概念，但下面的例子容易了解：
到這裡，清楚二者之前并無明顯聯系，我們回到任務的主線——
吉布斯現象（Gibbs phenomenon)
維基百科的描述是：
在工程應用時常用有限正弦項正弦波疊加逼近原周期信号。所用的諧波次數N的大小決定逼近原波形的程度，N增加，逼近的精度不斷改善。但是由于對于具有不連續點的周期信号會發生一種現象：當選取的傅裡葉級數的項數N增加時，合成的波形雖然更逼近原函數，但在不連續點附近會出現一個固定高度的過沖，N越大，過沖的最大值越靠近不連續點，但其峰值并不下降，而是大約等于原函數在不連續點處跳變值的9%，且在不連續點兩側呈現衰減振蕩的形式。
此前對現象的大緻含義了解，這次留意到的是關鍵詞：不連續點。
英文版的詞條頁面内容很多，顯然對這一現象的研究已非時日，我選取了最容易入手的一個方向：通過頁面提供的使用python代碼的優化，開啟記錄。
Fake Nodes
這是一種新的插值方法，消減吉布斯現象是其應用
插值定義(interpolation)
前兩天結束的美賽中隊友使用了這一關鍵詞描述模型，我初次聽聞。他的使用在于将資料集中以1 ∘ ^{\circ} ∘為分度值的溫度資料進行插值處理，以擷取連續位置的溫度值。
而這便是插值的的定義：
對于一組離散的資料點（稱為結點，是否正是與方法中的Nodes相對應？） ( x i , y i ) , i = 1 , 2 , . . . , n (x_{i},y_{i}),i=1,2,...,n (xi,yi),i=1,2,...,n, 對于 x ( x ≠ x i ) , i = 1 , 2 , . . . , n x(x\neq x_{i}),i=1,2,...,n x(x=xi),i=1,2,...,n ,求 x x x 對應的 y y y 值。
matlab代碼
%傳入離散的資料點
x=[0.0 0.1 0.195 0.3 0.401 0.5]
y=[0.39849 0.39695 0.39142 0.38138 0.36812 0.35206]
%确定a所對應的值T
a=.25;
%繪制離散資料點圖像
figure;
plot(x,y,'*-')
hold on
interpType='linear'
%不同的插值方法
T=interp1(x,y,a,interpType); %線性插值
plot(a,T,'go')
text(a-0.2,T,[interpType,':T','=',num2str(T)],'color','g')
hold on;
interpType='nearest'% 兩點插值
T=interp1(x,y,.25,interpType)
%(傳回結果T=0.3814)
plot(a,T,'ro')
text(a+0.1,T,[interpType,':T','=',num2str(T)],'color','r')
hold on;
interpType='spline'% 三次樣條插值
T=interp1(x,y,a,interpType)
%(傳回結果T =0.3867)
plot(a,T,'bo')
text(a+0.1,T,[interpType,':T','=',num2str(T)],'color','b')
hold on;
interpType='pchip'%三次插值
T=interp1(x,y,a,interpType)
plot(a,T,'ko')
text(a+0.1,T+0.003,[interpType,':T','=',num2str(T)],'color','k')
hold on;
title('插值函數demo')
結果：
對插值有了較為感性的認識。
龍格現象（Runge’s phenomenon）
運作了一遍關于Gibbs的notebook，内心毫無波瀾，來到前一個檔案Runge，且見開篇提到：
It is well known that the polynomial interpolation is unstable over equispaced samples, and stable the Chebyshev-Lobatto (CL) nodes.
衆所周知 多項式插值法應用于平均分布的樣本是不穩定的，而在Chebyshev-Lobatto (CL)類型結點上的應用是穩定的
OK,fine.我的确連所謂的衆所周知都一無所知，終于明白為何曾經高中英語作文時老師要求我們不要常用此句式，殺傷力【IIIIIIII】100%
繼續探索吧，學海無涯。
定義
龍格現象是對于平均分布的樣本應用多項式插值時，區間邊緣處出現震蕩問題。維基百科中确提到這與傅裡葉變換中的吉布斯現象是相似的，那有必要了解多項式的插值如何進行。
多項式插值(Polynomial interpolation)
前期提醒，其實應用僅僅涉及了線性代數的知識，耐心看完就會了解的。
首先定理如下（回顧之前的插值定義）：
Given a set of n + 1 n + 1 n+1 data points ( x i , y i ) (x_i, y_i) (xi,yi) where no two x i x_i xi are the same, one is looking for a polynomial p of degree at most n n n with the property
-
試譯：對于 n + 1 n + 1 n+1 個資料點 ( x i , y i ) (x_i, y_i) (xi,yi) ，其中 x i x_i xi 各異，尋找 n n n 階多項式 p p p 滿足：
p ( x i ) = y i , i = 0 , … , n . ( 1 ) {\displaystyle p(x_{i})=y_{i}, \qquad i=0,\ldots ,n.} \qquad(1) p(xi)=yi,i=0,…,n.(1)
根據唯一可解性，此多項式存在且唯一，可由範德蒙矩陣證明。
步驟
列出多項式 p ( x ) p(x) p(x):
p ( x ) = a n x n + a n − 1 x n − 1 + ⋯ + a 2 x 2 + a 1 x + a 0 . ( 2 ) p(x) = a_n x^n + a_{n-1} x^{n-1} + \cdots + a_2 x^2 + a_1 x + a_0. \qquad(2) p(x)=anxn+an−1xn−1+⋯+a2x2+a1x+a0.(2)
将一個資料點 ( x i , y i ) (x_i,y_i) (xi,yi) 代入式 ( 2 ) (2) (2)中得：
p ( x i ) = a n x i n + a n − 1 x i n − 1 + ⋯ + a 2 x i 2 + a 1 x i + a 0 . ( 3 ) p(x_i) = a_n x_i^n + a_{n-1} x_i^{n-1} + \cdots + a_2 x_i^2 + a_1 x_i + a_0. \qquad(3) p(xi)=anxin+an−1xin−1+⋯+a2xi2+a1xi+a0.(3)
由定理可得 ( 1 ) (1) (1)：
a n x i n + a n − 1 x i n − 1 + ⋯ + a 2 x i 2 + a 1 x i + a 0 = y i . ( 3 ) a_n x_i^n + a_{n-1} x_i^{n-1} + \cdots + a_2 x_i^2 + a_1 x_i + a_0=y_i. \qquad(3) anxin+an−1xin−1+⋯+a2xi2+a1xi+a0=yi.(3)
表示為矩陣形式是：
[ x i n x i n − 1 x i n − 2 … x i 1 ] [ a n a n − 1 ⋮ a 0 ] = [ y i ] . \begin{bmatrix} x_i^n & x_i^{n-1} & x_i^{n-2} & \ldots & x_i & 1 \end{bmatrix} \begin{bmatrix} a_n \\ a_{n-1} \\ \vdots \\ a_0 \end{bmatrix} = \begin{bmatrix} y_i \end{bmatrix}. [xinxin−1xin−2…xi1]⎣⎢⎢⎢⎡anan−1⋮a0⎦⎥⎥⎥⎤=[yi].
将每個資料點代入可得矩陣。
[ x 0 n x 0 n − 1 x 0 n − 2 … x 0 1 x 1 n x 1 n − 1 x 1 n − 2 … x 1 1 ⋮ ⋮ ⋮ ⋮ ⋮ x n n x n n − 1 x n n − 2 … x n 1 ] [ a n a n − 1 ⋮ a 0 ] = [ y 0 y 1 ⋮ y n ] . \begin{bmatrix} x_0^n & x_0^{n-1} & x_0^{n-2} & \ldots & x_0 & 1 \\ x_1^n & x_1^{n-1} & x_1^{n-2} & \ldots & x_1 & 1 \\ \vdots & \vdots & \vdots & & \vdots & \vdots \\ x_n^n & x_n^{n-1} & x_n^{n-2} & \ldots & x_n & 1 \end{bmatrix} \begin{bmatrix} a_n \\ a_{n-1} \\ \vdots \\ a_0 \end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_n \end{bmatrix}. ⎣⎢⎢⎢⎡x0nx1n⋮xnnx0n−1x1n−1⋮xnn−1x0n−2x1n−2⋮xnn−2………x0x1⋮xn11⋮1⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡anan−1⋮a0⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡y0y1⋮yn⎦⎥⎥⎥⎤.
學過線性代數後對這個形式就很熟悉了，當然和正常的符号表示有所差別，具體情況具體分析哦：
X A = Y XA=Y XA=Y
是以，求解多謝式 p p p 轉化為求解矩陣 A A A:
A = X − 1 Y A=X^{-1}Y A=X−1Y
實操
在日前觀看的Matlab教學視訊中（我有一個想法，用Matlab實作數理函數與特殊方程中的實體模型仿真，說法是什麼？挖坑？占坑？為啥要坑自己：（，了解到Matlab正是Matrix Lab的縮寫，用其處理矩陣乃上上之選。
資料點
資料點仍然選用上一節的範例：
x=[0.0 0.1 0.195 0.3 0.401 0.5]
y=[0.39849 0.39695 0.39142 0.38138 0.36812 0.35206]
矩陣 X X X
用循環的方法得到矩陣 X X X
- Matlab知識點：矩陣大小語句
size()
size(x)
的輸出結果是：
1 6 1 \qquad 6 16
那我們想得到第二個數值——矩陣的列數—— 6 6 6，則應輸入：
size(x,2)
-
Matlab知識點：編号起于 1 1 1
與常用的程式設計語言不同，這裡下标是從1開始
由于下标的問題，處理資料讓人稍有困惑。注意，這裡的 m m m代表 n + 1 n+1 n+1個資料點， n n n 即原公式中的 n n n（雞肋）
m=size(x,2)
n=m-1
X=zeros(m,m)
for i=1:m
for j=1:m
X(i,j)=x(i)^(m-j)
end
end
求解矩陣 X X X
- Matlab知識點：矩陣的逆
inv(X)
- Matlab知識點：矩陣的轉置
y'
%求多項式系數
a=inv(X)*y'
完整代碼及可視化展示
x=[0.0 0.1 0.195 0.3 0.401 0.5]
y=[0.39849 0.39695 0.39142 0.38138 0.36812 0.35206]
m=size(x,2)
n=m-1
X=zeros(m,m)
for i=1:m
for j=1:m
X(i,j)=x(i)^(m-j)
end
end
%求多項式系數
a=inv(X)*y'
d=0
for t=0:.001:0.5
d=d+1
%t=.25
sum=0
for i=1:m
sum=sum+a(i)*t^(m-i)
end
myinter(d)=sum
end
t=0:.001:0.5
plot (t,myinter)
hold on;
plot (x,y,'x-')
legend('myinter','data point')
總結
自己漸漸深入問題，拓展了知識面，實作效果不錯，正向激勵。寫完文章後作為第一階段，遙想這個問題需要多少階段……
心中無限滿足。