? ? ? ?在數(shù)值方法中，我們對時(shí)間、空間都做了有限寬度劃分，在求解問題時(shí)，我們得到物理量隨特定時(shí)刻、位置的分布（采樣是無法做到連續(xù)的）。于是，我們經(jīng)常會遇到這樣的問題：如何根據(jù)已知的采樣點(diǎn)數(shù)值，推算或者推測未知點(diǎn)的數(shù)值？

? ? ? ?當(dāng)然，如果我們有問題的解析解，比如物理量隨時(shí)間、空間變化的函數(shù)，直接代入即可。多數(shù)時(shí)候，我們采用數(shù)值方法是因?yàn)樵搯栴}不存在解析解或者說不存在簡單的函數(shù)形式，此時(shí)我們需要借助數(shù)據(jù)插值來解決。如下圖所示，我們知道體系能量隨晶格常數(shù)的變化（只有幾個(gè)點(diǎn)的數(shù)據(jù)），通過插值可以找到能量最低的位置；在比較能量和體積的關(guān)系時(shí)，通過插值可以找到兩種晶格轉(zhuǎn)變的臨界點(diǎn)。在二維問題中，插值可以提高圖片的分辨率，比如我們要掃描原子吸附的勢能面，可以計(jì)算少數(shù)位置的吸附能，得到比較全面的勢能面。

（1）??從線性插值到拉格朗日插值

? ? ? 最簡單的情況，如果知道平面上的兩個(gè)點(diǎn)(x_1,y_1)和(x_2，y_2)，我們可以假設(shè)方程形式為y=kx+b，解出k和b，最后寫成比較對稱的形式：當(dāng)x=x_1時(shí)，第一項(xiàng)是y_1，第二項(xiàng)是0，x=x_2時(shí)，第一項(xiàng)是0，第二項(xiàng)是y_2。這樣很容易推廣到N個(gè)點(diǎn)時(shí)對應(yīng)的N-1次多項(xiàng)式，也就是拉格朗日插值?？梢钥吹剑逯岛蟮那€通過之前已知的數(shù)據(jù)點(diǎn)，而未知的數(shù)據(jù)點(diǎn)可以根據(jù)多項(xiàng)式計(jì)算即可。注意，一般未知的自變量在已知點(diǎn)的“內(nèi)部”時(shí)結(jié)果要可靠，外推有風(fēng)險(xiǎn)。

? ? ? ?具體程序見附錄A和B。這里以A為例，x和y是已知的數(shù)據(jù)，我們用了cos函數(shù)的幾個(gè)點(diǎn)，一般是實(shí)驗(yàn)測量的數(shù)據(jù)。xx是自變量最小值到最大值之間的一個(gè)隨機(jī)數(shù)，根據(jù)公式，指標(biāo)i要原始數(shù)據(jù)所有點(diǎn)循環(huán)，一個(gè)有N個(gè)點(diǎn)。對于給定的i，j也要從1到N循環(huán)，但要把i去掉，所以我們用了setdiff這個(gè)集合運(yùn)算，再結(jié)合之前的prod函數(shù)計(jì)算連乘即可，L的第i項(xiàng)是其中一個(gè)分量，最后結(jié)果是sum(L) 求和。對于程序B，xx是已知數(shù)據(jù)自變量最小值到最大值的區(qū)間，有100個(gè)點(diǎn)，對xx循環(huán)，然后L(i,k)就是給定xx(k)插值時(shí)的某一分量，最后可以畫出L(i,:)隨xx的變化，L的某一行對應(yīng)的函數(shù)只通過原始數(shù)據(jù)的某一個(gè)點(diǎn)。

（2）Octave的interp1、spline函數(shù)

? ? ??在平時(shí)的使用中，我們掌握插值函數(shù)的語法即可，這里主要介紹Octave里面的1維、2維插值，分別對應(yīng)interp1 和 interp2 函數(shù)。在下面程序（附錄C）中，a的兩列分別是已知數(shù)據(jù)的自變量、因變量，如果我們直接用直線連接就是藍(lán)線。x=linspace(min(a(:,1)),max(a(:,1)),200);?是把自變量區(qū)間生成200個(gè)點(diǎn)，y=interp1(a(:,1),a(:,2),x,'spline'); 是調(diào)用插值函數(shù) interp1，注意語法，里面有4個(gè)參數(shù)，分別對應(yīng)已知數(shù)據(jù)的自變量、因變量、我們關(guān)心的自變量、插值方法。這里的'spline'是樣條插值方法，對應(yīng)分段多項(xiàng)式，并有一定光滑的連接要求，具體思路可以查閱相應(yīng)書籍。其中spline在Octave中也是獨(dú)立的插值函數(shù)，第9行可以換成y=spline(a(:,1),a(:,2),x);得到的結(jié)果是相同的。y是插值后的因變量（紅線），通過min函數(shù)查找最小值，[q,qq]=min(y)，這里q返回最小值，qq返回最小值所在位置對應(yīng)的整數(shù)，所以x(qq)就是讓y最小的自變量。

? ? 在附錄D中，我們有兩組曲線，分別保存在a和b矩陣，可以分別插值。為了求交點(diǎn)，我們需要用相同的自變量區(qū)間，這里 x=linspace(max([min(a(:,1)) min(b(:,1))]),min([max(a(:,1)) max(b(:,1))])); 取兩組數(shù)據(jù)的共用部分，保證數(shù)據(jù)在內(nèi)部，避免外推的誤差。用ay=spline(a(:,1),a(:,2),x); by=spline(b(:,1),b(:,2),x);分別插值后，計(jì)算ay和by差的絕對值，因?yàn)榻稽c(diǎn)對應(yīng)的是ay和by幾乎相同。[q,qq]=min(abs(ay-by))；根據(jù)x(qq) 得到交點(diǎn)的自變量。

（3）二維插值和圖像處理

? ? ?Octave的二維插值函數(shù)是interp2，但要注意自變量是矩陣，因?yàn)榇藭r(shí)是矩陣之間的映射。在附錄E中，x和y是1維向量，通過for循環(huán)計(jì)算得到 V(i,j)=exp(-((i-6)^2+(j-6)^2)/10); 因?yàn)榫W(wǎng)格稀疏，用image得到圖片有明顯的方格。先用xi=1:0.02:N; yi=1:0.02:N; 對x、y加密（增加劃分點(diǎn)），然后[xi,yi]=meshgrid(xi,yi); 得到對應(yīng)加密后的網(wǎng)格，最后通過zi=interp2(x,y,temps,xi,yi,'cubic'); 實(shí)現(xiàn)二維插值，這里'cubic'也是可以選擇不同方法，help interp2查看說明即可。

? ? ? 附錄F是對圖片的處理，這里選擇了卡通圖片，先把采樣點(diǎn)變少，然后看插值后的效果，能實(shí)現(xiàn)模糊的平滑，但是無法變清晰，這是因?yàn)樗惴ㄊ腔诠饣瘮?shù)，沒有“銳化”的部分，將來有機(jī)會再用神經(jīng)網(wǎng)絡(luò)和多圖片訓(xùn)練嘗試一下。

程序附錄：

A. 拉格朗日方法--單點(diǎn)

clear

x=linspace(0,2*pi,8);

y=cos(x);

N=length(x);

plot(x,y,'o')

xx=min(x)+rand*(max(x)-min(x));

for i=1:N

? j=setdiff(1:N,i);

? L(i)=y(i)*prod(xx-x(j))/prod(x(i)-x(j));

end

sum(L)

B. 拉格朗日方法--區(qū)間

clear

x=linspace(0,2*pi,6);

y=cos(x);

N=length(x);

xx=linspace(min(x),max(x));

for k=1:length(xx)

? for i=1:N

? j=setdiff(1:N,i);

? L(i,k)=y(i)*prod(xx(k)-x(j))/prod(x(i)-x(j));

end

end

plot(x,y,'ob')

hold on?

plot(xx,sum(L),'r')

figure

plot(x,y,'ob')

hold on?

for i=1:N

? plot(xx,L(i,:))

end

C. 1維插值求最低點(diǎn)

clear

a=[26.8892? ?-3.2733

? ?36.8850? ?-3.8328

? ?49.0939? ?-3.5725

? ?63.7373? ?-3.0755];

plot(a(:,1),a(:,2),'o-')

hold on?

x=linspace(min(a(:,1)),max(a(:,1)),200);

y=interp1(a(:,1),a(:,2),x,'spline');

plot(x,y,'r')??

[q,qq]=min(y)

x(qq)

?

D.? 1維插值求曲線交點(diǎn)

clear

a=[26.8892? ?-3.2733

? ?36.8850? ?-3.8328

? ?49.0939? ?-3.5725

? ?63.7373? ?-3.0755];

plot(a(:,1),a(:,2),'ob')

ax=linspace(min(a(:,1)),max(a(:,1)),200);

ay=spline(a(:,1),a(:,2),ax);

hold on?

plot(ax,ay,'b')

b=[? ?20.2589? ?-3.1978

? ?27.7900? ?-3.7763

? ?36.9885? ?-3.5219

? ?48.0211? ?-3.0122];

plot(b(:,1),b(:,2),'or')

bx=linspace(min(b(:,1)),max(b(:,1)),200);

by=spline(b(:,1),b(:,2),bx);

hold on?

plot(bx,by,'r')

x=linspace(max([min(a(:,1)) min(b(:,1))]),min([max(a(:,1)) max(b(:,1))]));

ay=spline(a(:,1),a(:,2),x);

by=spline(b(:,1),b(:,2),x);

[q,qq]=min(abs(ay-by))；

x(qq)

E. 二維插值舉例

clear

N=11;x=1:N;y=1:N;

for i=1:length(x)

for j=1:length(y)

V(i,j)=exp(-((i-6)^2+(j-6)^2)/10);

end

end

temps=V;

figure(1)

image(V*50)

colormap(jet)

figure(2);

xi=1:0.02:N;

yi=1:0.02:N;

[xi,yi]=meshgrid(xi,yi);

zi=interp2(x,y,temps,xi,yi,'cubic');

image(zi*50)

colormap(jet)

F. 圖片插值處理

clear

a=imread('a.jpg');

subplot(1,3,1)

imshow(a)

a=a(1:20:size(a,1),1:20:size(a,2),:);

subplot(1,3,2)

imshow(a)

xi=1:0.1:size(a,2);

yi=1:0.1:size(a,1);

[xi,yi]=meshgrid(xi,yi);

zi(:,:,1)=interp2(1:size(a,2),1:size(a,1),a(:,:,1),xi,yi,'cubic');

zi(:,:,2)=interp2(1:size(a,2),1:size(a,1),a(:,:,2),xi,yi,'cubic');

zi(:,:,3)=interp2(1:size(a,2),1:size(a,1),a(:,:,3),xi,yi,'cubic');

subplot(1,3,3)

imshow(zi)