埃尔米特插值c语言_计算物理基础：第七章-插值与拟合

参考北京师范大学的《计算物理基础》

第七章-插值与拟合

计算物理基础_中国大学MOOC(慕课)​www.icourse163.org

1.拉格朗日插值法

已知函数f(x)的某些离散值，希望用函数g来逼近f，就会需要补充一些点的值，使其形成光滑的曲线，这种方法叫做插值法。

1.1插值分类

• 按插值范围：内插、外插
• 按插值方法：多项式插值、样条插值、曲线插值、分段插值

1.2拉格朗日插值法

• 一次插值函数

• 二次插值函数

• 一般情况

给定函数

在

个节点的函数值，要求出一个不超过n次的插值多项式。

代码实现：

function test
x=1:6;
y=[16,18,21,17,15,12];
u=0.75:0.05:6.25;
v=polyinterp(x,y,u)
plot(x,y,'*:',u,v,'o')  %*：原数据点，o：新插入点，o：插值曲线，虚线：分段线性插值

function v=polyinterp(x,y,u) %x、y0是原有数据点，u是新插入的点
n=length(x); 
v=zeros(size(u));
for k=1:n 
    w=ones(size(u));
    for j=[1:k-1 k+1:n]
        w=(u-x(j))./(x(k)-x(j)).*w;
    end
v=v+w*y(k);
end

2.插值指令的用法

2.1拉格朗日插值法与分段线性插值法的优缺点比较：

它们是插值方法的两个极端，需要改进：

2.2改进思路：

• 分段插值：将数据分成小区域插值；
• 埃尔米特插值：使分段曲线光滑连接，连接点有一阶（二阶）导数。

不同方法对导数要求不同，指令pchip（保形分段三次插值：一阶导数连续）和spline（样条插值：二阶导数连续）是两种方法。

如何实现连接分段的插值曲线光滑且具有连续导数呢？

2.3 分段三次埃尔米特插值

指令对端点导数值选择

matlab插值指令用法：x1是新插入的值

• 实例：以正弦函数的10个值作为原始数据（红点），用4种方法（最近邻插值、线性插值、分段三次样条插值、保形分段三次插值）算41个插值点并画线。

x=0:10;
y=sin(x);
xi=0:.15:10;%插值点
yi1=interp1(x,y,xi,'nearest');
yi2=interp1(x,y,xi,'linear');
yi3=interp1(x,y,xi,'spline');
yi4=interp1(x,y,xi,'pchip');

subplot(2,2,1)
plot(x,y,'r*',xi,yi1)
subplot(2,2,2)
plot(x,y,'r*',xi,yi2)
subplot(2,2,3)
plot(x,y,'r*',xi,yi3)
subplot(2,2,4)
plot(x,y,'r*',xi,yi4)

可见，样条插值的曲线最光滑。

3.曲线拟合

3.1曲线拟合

曲线拟合：也叫函数逼近，是另一种插值方法。

本节：用最小二乘法作曲线拟合。

3.2曲线拟合与拉格朗日插值比较：

• 拉格朗日插值：要求函数通过数据点、函数与数据点偏差为零；
• 曲线拟合：不要求函数通过数据点、可用存在偏差，反映逼近程度。

3.3 插值函数与数据点的三种偏差准则

那么在曲线拟合中，如何估算插值函数与数据点的偏差呢？有以下三种偏差判断准则：（第三种是使用最广泛、计算比较简单的方法）

3.4 用多项式作为拟合函数所采用的计算公式：

• 设所求插值函数是一个m次多项式
  

• 选取多项式的系数
  使得在
  点有

最小，由多元函数取极值的条件得：

，将其整理，得

多项式拟合指令polyfit：matlab将多项式拟合开发成专用的指令，在指令的程序中使用了还能先进的算法来求解正则方程：

ployfit(X,Y,N):最小二乘法多项式拟合
polyval(P,X):计算多项式的值

X、Y：数据点及其函数值
N：拟合多项式的最高幂次
P：被求值的多项式系数，按降幂排列

例1：

t=0:0.5:5;
s=1*t+0.2.*t.*t;%给定二次曲线数据点
s1=[0.5 -0.18 -0.01 0.13 0.1...%随机数列
0.31 -0.22 -0.31 0.2 0.4 -0.14];
ss=s+s1;%设定误差值
P=polyfit(t,ss,2) %二次曲线拟合
y=polyval(P,t);%计算拟合函数的值
plot(t,y,t,ss,'r*')

%输出结果
P =

    0.2261    0.8433    0.2345

拟合曲线没有通过每一点，但是得到的的曲线很光滑。

例2：对比推导的正则方程组的解与用指令求得的解

————————方程组解————————

设拟合多项式为

正则方程组为

左除法求得方程的解为

代码：

A=[9,0,3.75;0,3.75,0;3.75,0,2.7656]; 
b=[18.1732;8.4842;7.6173];
X=Ab

%输出结果
X =

    2.0036
    2.2625
    0.0375

————————指令求解————————

用指令求解：

x=[-1 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1];
y=[-0.2209 0.3295 0.8826 1.4392 2.0003 2.5645 3.1334 3.7601 4.2836]; 
yi=polyfit(x,y,2)


%输出结果
yi =

    0.0379    2.2624    2.0034

可见，用指令得到的系数与用方程组解求得的系数相近。

4.曲线拟合工具箱

matlab中对曲线拟合采用的是工具箱方法来展现各种指令、进行拟合。

功能简介：

操作步骤：

例1：

首先在实时脚本文件处输入数据t和I

>> t=0.2:0.1:0.8

t =

    0.2000    0.3000    0.4000    0.5000    0.6000    0.7000    0.8000

>> I=[3.16,2.38,1.75,1.34,1.00,0.74,0.56]

I =

    3.1600    2.3800    1.7500    1.3400    1.0000    0.7400    0.5600

输入cftool指令，启动曲线拟合工具箱，其页面布局有：数据输入、函数选择、拟合结果的偏差和公式，输出图形，将x data置为t、y data置为I，函数选择“Exponential”（指数函数），即可输出函数方程、图像、a和b的值，如下图所示

例2：

首先在实时脚本文件处输入数据x和y

>> x=[40:10:400]/180*pi;%弧度表示
>> y=[0.02 0.04 0.1 0.21 0.39 0.46 0.63 0.72 0.81 0.91 0.96 0.92 0.84 0.72 0.5 0.33 0.18 0.08...
0.03 0.03 0.09 0.22 0.39 0.57 0.72 0.88 0.99 1.06 0.98 0.83 0.67 0.53 0.41 0.28 0.17 0.06 0.03];

将x data置为x、y data置为y，函数选择“Custom Equation”（自定义函数），在填框中输入函数形式“a*(cos(b*x+c))^2”，即可输出函数方程自变量a、b、c的值及拟合曲线，如下图所示

注：也可通过自动拟合“Auto fit”进行拟合、多次拟合。从而得到曲线和公式中的参数。

2020.12.12