使用背景:在工程实践和科学实验中, 常常需要从一组观测数据之中找到自变量与因变量间的函数关系,一般可用一个近似函数来表示。函数的产生办法因观测数据和要求的不同而各异,数据插值与拟合是两种常用的方法。
函数插值来源于函数的以下问题:
只知道函数在某区间有定义且已得到区间内的一些离散的点的值,希望能用简单的表达式近似给出函数在此区间上的整体描述,并能与已知离散点上的值相等。
language插值代码:
function L=language(x,y,x0)
n=length(x);
syms t;
L=0;
for k=1:n
lx=1;
for j=1:k-1
lx=lx*(t-x(j))/(x(k)-x(j));
end
for j=k+1:n
lx=lx*(t-x(j))/(x(k)-x(j));
end
L=L+y(k)*lx;
end
if nargin==3
L=subs(L,'t',x0);
else
L=collect(L);
L=vpa(L,6);
end
牛顿插值法:
function N=newton(x,y,x0)
n=length(x);
syms t;
N=y(1);
f=0;
lx=1;
for i=1:n-1
for j=i+1:n
f(j)=(y(j-1)-y(j))/(x(j-i)-x(j));
end
lx=lx*(t-x(i));
N=N+ f(i+1)*lx;
y=f;
end
if nargin==3
N=subs(N,'t',x0);
else
N=collect(N);
N=vpa(N,6);
end
函数逼近与曲线拟合:
函数逼近:
在区间[a, b]上已知一连续函数 f (x), 如果 f (x)的表达式太过复杂, 可用一简单函数去近似 f (x), 这就是函数逼近问题。
曲线拟合:
如果 f (x)的表达式未知, 只知道描述 f (x)的一条曲线, 这就是曲线拟合问题。
和插值问题不同的是, 逼近和拟合并不要求逼近函数在已知点上的值一定得等于原函数的函数值, 而是按照某种标准使得两者的差值达到最小。
切比雪夫逼近:
当一个连续函数定义在区间[-1, 1]上时, 它可以展开成切比雪夫级数叠加形式:
function f=chebyshev(fun,k,x0)
syms t;
T(1:k+1)=t;
T(1)=1; T(2)=t;
c(1:k+1)=sym(0);
fun=subs(fun,findsym(sym(fun)),t);
c(1)=1*int(fun*T(1)/sqrt(1-t^2),t,-1,1)/pi;
c(2)=2*int(fun*T(2)/sqrt(1-t^2),t,-1,1)/pi;
f=c(1)*T(1)+c(2)*T(2);
for n=3:k+1
T(n)=2*t*T(n-1)-T(n-2);
c(n)=2*int(fun*T(n)/sqrt(1-t^2),t,-1,1)/pi;
f=f+c(n)*T(n);
end
if nargin==3
f=subs(f,'t',x0); f=vpa(f,6);
else
f=collect(f); f=vpa(f,6);
end
勒让德逼近:
勒让德逼近也要求被逼近函数定义在区间[-1, 1]上, 且被逼近函数可以展开成叠加形式:
function f=legendre(fun,k,x0)
syms t;
P(1:k+1)=t;
P(1)=1;P(2)=t;
c(1:k+1)=sym(0);
fun=subs(fun,findsym(sym(fun)),t);
c(1)=1*int(fun*P(1),t,-1,1)/2;
c(2)=3*int(fun*P(2),t,-1,1)/2;
f=c(1)*P(1)+c(2)*P(2);
for n=3:k+1
P(n)=((2*n-3)*t*P(n-1)-(n-2)*P(n-2))/(n-1);
c(n)=(2*n-1)*int(fun*P(n),t,-1,1)/2;
f=f+c(n)*P(n);
end
if nargin==3
f=subs(f,'t',x0); f=vpa(f,6);
else
f=collect(f); f=vpa(f,6);
end
傅里叶逼近:
在很多情形下, 我们所遇到的函数 f (x)只在某个区间上才有意义, 例如固定在 x= -l 和 l 处的弦振动问题, 弦上个点的振动只在区间[-l, l]上才有意义, 此时我们可以把它展开为傅立叶级数叠加形式:
function [a0,a,b]=fourier(fun,T,k)
syms t;
fun=subs(fun,findsym(sym(fun)),t);
a0=int(fun,t,-T/2,T/2)/T;
for n=1:k
a(n)=int(fun*cos(2*n*pi*t/T),t,-T/2,T/2)*2/T;
b(n)=int(fun*sin(2*n*pi*t/T),t,-T/2,T/2)*2/T;
end
a=vpa(a,6);
b=vpa(b,6);
多项式拟合:
对于给定的实验数据点(xi, yi) (i=1, 2, …, N), 可构造m次多项式
function a=multifit(x,y,m)
N=length(x);
c(1:2*m+1)=0;
b(1:m+1)=0;
for k=1:2*m+1
for i=1:N
c(k)=c(k)+x(i)^(k-1);
if k<=m+1
b(k)=b(k)+y(i)*x(i)^(k-1);
end
end
end
for s=1:m+1
C(s,:)=c(s:m+s);
end
a=C^(-1)*b';
最小二乘法拟合:
最小二乘法拟合也是科学实验中经常使用的方法, 它的具体操作过程是从一组试验数据(xi, yi)中拟合出函数关系 y=f (x), 拟合的标准是使得 ( f(xi) - yi)的平方取极小值。
命令: polyfit(x,y,n)