计算物理基础-数据插值
在数值方法中，我们对时间、空间都做了有限宽度划分，在求解问题时，我们得到物理量随特定时刻、位置的分布（采样是无法做到连续的）。于是，我们经常会遇到这样的问题：如何根据已知的采样点数值，推算或者推测未知点的数值？
当然，如果我们有问题的解析解，比如物理量随时间、空间变化的函数，直接代入即可。多数时候，我们采用数值方法是因为该问题不存在解析解或者说不存在简单的函数形式，此时我们需要借助数据插值来解决。如下图所示，我们知道体系能量随晶格常数的变化（只有几个点的数据），通过插值可以找到能量最低的位置；在比较能量和体积的关系时，通过插值可以找到两种晶格转变的临界点。在二维问题中，插值可以提高图片的分辨率，比如我们要扫描原子吸附的势能面，可以计算少数位置的吸附能，得到比较全面的势能面。
（1） 从线性插值到拉格朗日插值
最简单的情况，如果知道平面上的两个点(x_1,y_1)和(x_2，y_2)，我们可以假设方程形式为y=kx+b，解出k和b，最后写成比较对称的形式：当x=x_1时，第一项是y_1，第二项是0，x=x_2时，第一项是0，第二项是y_2。这样很容易推广到N个点时对应的N-1次多项式，也就是拉格朗日插值。可以看到，插值后的曲线通过之前已知的数据点，而未知的数据点可以根据多项式计算即可。注意，一般未知的自变量在已知点的“内部”时结果要可靠，外推有风险。
具体程序见附录A和B。这里以A为例，x和y是已知的数据，我们用了cos函数的几个点，一般是实验测量的数据。xx是自变量最小值到最大值之间的一个随机数，根据公式，指标i要原始数据所有点循环，一个有N个点。对于给定的i，j也要从1到N循环，但要把i去掉，所以我们用了setdiff这个集合运算，再结合之前的prod函数计算连乘即可，L的第i项是其中一个分量，最后结果是sum(L) 求和。对于程序B，xx是已知数据自变量最小值到最大值的区间，有100个点，对xx循环，然后L(i,k)就是给定xx(k)插值时的某一分量，最后可以画出L(i,:)随xx的变化，L的某一行对应的函数只通过原始数据的某一个点。
（2）Octave的interp1、spline函数
在平时的使用中，我们掌握插值函数的语法即可，这里主要介绍Octave里面的1维、2维插值，分别对应interp1 和 interp2 函数。在下面程序（附录C）中，a的两列分别是已知数据的自变量、因变量，如果我们直接用直线连接就是蓝线。x=linspace(min(a(:,1)),max(a(:,1)),200); 是把自变量区间生成200个点，y=interp1(a(:,1),a(:,2),x,'spline'); 是调用插值函数 interp1，注意语法，里面有4个参数，分别对应已知数据的自变量、因变量、我们关心的自变量、插值方法。这里的'spline'是样条插值方法，对应分段多项式，并有一定光滑的连接要求，具体思路可以查阅相应书籍。其中spline在Octave中也是独立的插值函数，第9行可以换成y=spline(a(:,1),a(:,2),x);得到的结果是相同的。y是插值后的因变量（红线），通过min函数查找最小值，[q,qq]=min(y)，这里q返回最小值，qq返回最小值所在位置对应的整数，所以x(qq)就是让y最小的自变量。
在附录D中，我们有两组曲线，分别保存在a和b矩阵，可以分别插值。为了求交点，我们需要用相同的自变量区间，这里 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);分别插值后，计算ay和by差的绝对值，因为交点对应的是ay和by几乎相同。[q,qq]=min(abs(ay-by))；根据x(qq) 得到交点的自变量。
（3）二维插值和图像处理
Octave的二维插值函数是interp2，但要注意自变量是矩阵，因为此时是矩阵之间的映射。在附录E中，x和y是1维向量，通过for循环计算得到 V(i,j)=exp(-((i-6)^2+(j-6)^2)/10); 因为网格稀疏，用image得到图片有明显的方格。先用xi=1:0.02:N; yi=1:0.02:N; 对x、y加密（增加划分点），然后[xi,yi]=meshgrid(xi,yi); 得到对应加密后的网格，最后通过zi=interp2(x,y,temps,xi,yi,'cubic'); 实现二维插值，这里'cubic'也是可以选择不同方法，help interp2查看说明即可。
附录F是对图片的处理，这里选择了卡通图片，先把采样点变少，然后看插值后的效果，能实现模糊的平滑，但是无法变清晰，这是因为算法是基于光滑函数，没有“锐化”的部分，将来有机会再用神经网络和多图片训练尝试一下。
程序附录：
A. 拉格朗日方法--单点
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. 拉格朗日方法--区间
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维插值求最低点
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维插值求曲线交点
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)
