matlab代码--函数逼近

matlab代码--函数逼近
函数逼近
1.Chebyshev用切比雪夫多项式逼近已知函数
function f=Chebyshev(y,k,x0)
syms t;
T(1:k+1)=t;
T(1)=1;
T(2)=t;
c(1:k+1)=0.0;
c(1)=int(subs(y,findsym(sym(y)),sym('t'))*T(1)/sqrt(1-t^2),t,-1,1)/pi;
c(2)=2*int(subs(y,findsym(sym(y)),sym('t'))*T(2)/sqrt(1-
t^2),t,-1,1)/pi;
f=c(1)+c(2)*t;
for i=3:k+1
T(i)=2*t*T(i-1)-T(i-2);
c(i)=2*int(subs(y,findsym(sym(y)),sym('t'))*T(i)/sqrt(1-t^2),t,-1,1)/2;
f=f+c(i)*T(i);
f=vpa(f,6);
if(i==k+1)
if(nargin==3)
f=subs(f,'t',x0);
else
f=vpa(f,6);
end
end
end
2.Legendre用勒让德多项式逼近已知函数
function f=Legendre(y,k,x0)
syms t;
P(1:k+1)=t;
P(1)=1;
P(2)=t;
c(1:k+1)=0.0;
c(1)=int(subs(y,findsym(sym(y)),sym('t'))*P(1),t,-1,1)/2; c(2)=int(subs(y,findsym(sym(y)),sym('t'))*P(2),t,-1,1)/2; f=c(1)+c(2)*t;
for i=3:k+1
P(i)=((2*i-3)*P(i-1)*t-(i-2)*P(i-2))/(i-1);
c(i)=int(subs(y,findsym(sym(y)),t)*P(i),t,-1,1)/2;
f=f+c(i)*P(i);
if(i==k+1)
if(nargin==3)
f=subs(f,'t',x0);
else
f=vpa(f,6);
end
end
end
3.Pade用帕德形式的有理分式逼近已知函数
function f=Pade(y,n,x0)
syms t;
A=zeros(n,n);
q=zeros(n,1);
p=zeros(n+1,1);
b=zeros(n,1);
yy=0;
a(1:2*n)=0.0;
for(i=1:2*n)
yy=diff(sym(y),findsym(sym(y)),n);
a(i)=subs(sym(yy),findsym(sym(yy)),0.0)/factorial(i); end; for(i=1:n)
for(j=1:n)
A(i,j)=a(i+j-1);
end;
b(i,1)=-a(n+i);
end;
q=A\b;
p(1)=subs(sym(y),findsym(sym(y)),0.0);
for(i=1:n)
p(i+1)=a(n)+q(i)*subs(sym(y),findsym(sym(y)),0.0);
for(j=2:i-1)
p(i+1)=p(i+1)+q(j)*a(i-j);
end
end
f_1=0;
f_2=1;
for(i=1:n+1)
f_1=f_1+p(i)*(t^(i-1));
end
for(i=1:n)
f_2=f_2+q(i)*(t^i);
end
if(nargin==3)
f=f_1/f_2;
f=subs(f,'t',x0);
else
f=f_1/f_2;
f=vpa(f,6);
end
4.lmz用列梅兹算法确定函数的最佳一致逼近多项式function[coff,err]=lmz(func,m,a,b,eps)
if(nargin==4)
eps=1.0e-6;
end
syms v;
maxv=0.0;
max_x=a;%记录abs(f(x)-p(x))取最大值的x
for k=0:m
px(k+1)=power(v,k);
end%p（x）多项式
for i=1:m+2
x(i)=0.5*(a+b+(b-a)*cos(3.14159265*(m+2-i)/(m+1)));
fx(i)=subs(sym(func),findsym(sym(func)),x(i));
end%初始的x和f(x)
A=zeros(m+2,m+2);
for i=1:m+2
for j=1:m+1
A(i,j)=power(x(i),j-1);
end
A(i,m+2)=(-1)^i;
end
c=A\transpose(fx);%p(x)的初始系数
u=c(m+2);%算法中的u
tol=1;%精度
while(tol>eps)
t=a;
while(t<b)%此循环找出abs(f(x)-p(x))取最大值的x< p="">
t=t+0.05*(b-a)/m;
px1=subs(px,'v',t);
pt=px1*c(1:m+1);
ft=subs(sym(func),findsym(sym(func)),t);
if abs(ft-pt)>maxv
maxv=abs(ft-pt);
max_x=t;
end
end
if max_x>b
max_x=b;
end
%以下可参考算法的三个确定新点集的情况
if(a<=max_x)&&(max_x<=x(2))%第一种情况f0=subs(sym(func),findsym(sym(func)),x(2));
px1=subs(px,'v',x(2));
pt=px1*c(1:m+1);
d1=f0-pt;
fm=subs(sym(func),findsym(sym(func)),max_x);
pm1=subs(px,'v',max_x);
pm=pm1*c(1:m+1);
d2=fm-pm;
if d1*d2>0
x(2)=max_x;
end
else
if(x(m+1)<=max_x)&&(max_x<=b)%第二种情况
f0=subs(sym(func),findsym(sym(func)),x(m+1));
px1=subs(px,'v',x(m+1));
pt=px1*c(1:m+1);
d1=f0-pt;
fm=subs(sym(func),findsym(sym(func)),max_x); pm1=subs(px,'v',max_x);
pm=pm1*c(1:m+1);
d2=fm-pm;
if d1*d2>0
x(m+1)=max_x;
end
else%第三种情况
for i=2:m
if(x(i)<=max_x)&&(x(i+1)>=max_x)
index_x=i;
break;
end
end%找到max_x所在区间
f0=subs(sym(func),findsym(sym(func)),x(index_x)); px1=subs(px,'v',x(index_x));
pt=px1*c(1:m+1);
d1=f0-pt;
fm=subs(sym(func),findsym(sym(func)),max_x); pm1=subs(px,'v',max_x);
pm=pm1*c(1:m+1);
d2=fm-pm;
if d1*d2>0
x(index_x)=max_x;
end
end
end
for i=1:m+2%重新计算f（x）
fx(i)=subs(sym(func),findsym(sym(func)),x(i));
end
for i=1:m+2
for j=1:m+1
A(i,j)=power(x(i),j-1);
end
A(i,m+2)=(-1)^i;
end
c=A\transpose(fx);%重新计算p（x）的系数
tol=abs(c(m+2)-u);
u=c(m+2);
end
coff=c(1:m+1);
err=u;
5.ZJPF求已知函数的最佳平方逼近多项式
function coff=ZJPF(func,n,a,b)
C=zeros(n+1,n+1);
var=findsym(sym(func));
func=func/var;
for i=1:n+1
C(1,i)=(power(b,i)-power(a,i))/i;%算法中的C矩阵的第一行func=func*var;
d(i,1)=int(sym(func),var,a,b);%算法中的D向量的第一行end for i=2:n+1
C(i,1:n)=C(i-1,2:n+1);
f1=power(b,n+i);
f2=power(a,n+i);
C(i,n+1)=(f1-f2)/(n+i);%形成C矩阵
end
coff=C\d;%求解逼近多项式的系数6.FZZ用傅立叶级数逼近已知的连续周期函数
function[A0,A,B]=FZZ(func,T,n)
syms t;
func=subs(sym(func),findsym(sym(func)),sym('t')); A0=int(sym(func),t,-T/2,T/2)/T;
for(k=1:n)
A(k)=int(func*cos(2*pi*k*t/T),t,-T/2,T/2)*2/T;
A(k)=vpa(A(k),4);
B(k)=int(func*sin(2*pi*k*t/T),t,-T/2,T/2)*2/T;
B(k)=vpa(B(k),4);
end
7.DFF离散周期数据点的傅立叶逼近
function c=DFF(f,N)
c(1:N)=0;
for(m=1:N)
for(n=1:N)
c(m)=c(m)+f(n)*exp(-i*m*n*2*pi/N);
end
c(m)=c(m)/N;
end
8.SmartBJ用自适应分段线性法逼近已知函数function[node,err]=SmartBJ(func,a,b,maxtol) format long;
node(1)=a;
node(2)=b;
num=2;
if(b-a)<10
n=5;
else
n=10;
end
err=0;
bSign=1;
while(bSign)
bSign=0;
knode=node;
tnum=num;
insert_num=0;
for i=1:(tnum-1)
[mx,mv]=FindMX(func,knode(i),knode(i+1),n);
%找到区间[knode(i),knode(i+1)]上的误差最大的点
if mv>maxtol
%如果误差超过给定精度，在此点将区间[knode(i),knode(i+1)]分为两段d(1:(i+insert_num))=node(1:(i+insert_num));
d(i+insert_num+1)=mx;
num=num+1;
d((i+insert_num+2):num)=node((i+insert_num+1):(num-1));
node=d;
bSign=1;
insert_num=insert_num+1;
else
if(mv>err)
err=mv;%记录所有分段线性插值区间上的最大误差end
end
end
end
format short;
function[max_x,max_v]=FindMX(func,a,b,n)
format long;
eps=1.e-3;
max_v=0;
max_x=a;
fa=subs(sym(func),findsym(sym(func)),a);%左端点函数值fb=subs(sym(func),findsym(sym(func)),b);%右端点函数值step=n/5;
tol=1;
tmp=0;
while tol>eps
t=a;
for j=0:(n/step)%此循环找出取最大值的x
t=a+j*step*(b-a)/n;
pt=fa+(t-a)*(fb-fa)/(b-a);%线性插值得出的函数值
ft=subs(sym(func),findsym(sym(func)),t);
if abs(ft-pt)>max_v%abs(f(x)-p(x))
max_v=abs(ft-pt);%记录最大误差
max_x=t;%记录此点坐标
end
end
tol=abs(max_x-tmp);
tmp=max_x;
step=step/2;
end
format short;
9.SmartYTBJ用自适应样条逼近（第一类）已知函数function[node,err]=SmartYTBJ(func,a,b,y_s,y_e,maxtol) format long;
node(1)=a;
node(2)=b;
num=2;
if(b-a)<10
n=5;
else
n=10;
end
err=0;
bSign=1;
while(bSign)
bSign=0;
for l=1:num
y(l)=subs(subs(sym(func),findsym(sym(func)),node(l)));
end
knode=node;
tnum=num;
insert_num=0;
for i=1:(tnum-1)
pfx=ThrSample1(knode,y,y_s,y_e,((knode(i)+knode(i+1))/2));
%上式给出每个分段区间上的样条插值函数
[mx,mv]=FindMX(func,pfx,knode(i),knode(i+1),n);
%找到区间[knode(i),knode(i+1)]上的误差最大的点
if mv>maxtol
%如果误差超过给定精度，在此点将区间[knode(i),knode(i+1)]分为两段，即将此点加
%入结点数组中
d(1:(i+insert_num))=node(1:(i+insert_num));
d(i+insert_num+1)=mx;
num=num+1;
d((i+insert_num+2):num)=node((i+insert_num+1):(num-1));
node=d;
bSign=1;
insert_num=insert_num+1;
else
if(mv>err)
err=mv;%记录所有样条插值区间上的最大误差
end
end
end
end
format short;
function[max_x,max_v]=FindMX(func,pfx,a,b,n)
format long;
eps=1.e-3;
max_v=0;
max_x=a;
fa=subs(sym(func),findsym(sym(func)),a);%左端点函数值
fb=subs(sym(func),findsym(sym(func)),b);%右端点函数值
step=n/5;
tol=1;
tmp=0;
while tol>eps
t=a;
for j=0:(n/step)%此循环找出取最大值的x
t=a+j*step*(b-a)/n;
pt=subs(sym(pfx),findsym(sym(pfx)),t);%样条插值得出的函数值
ft=subs(sym(func),findsym(sym(func)),t);
if abs(ft-pt)>max_v%abs(f(x)-p(x))
max_v=abs(ft-pt);%记录最大误差
max_x=t;%记录此点坐标
end
end
tol=abs(max_x-tmp);
tmp=max_x;
step=step/2;
end
format short;
10.multifit离散试验数据点的多项式曲线拟合
function A=multifit(X,Y,m)
%A--输出的拟合多项式的系数
N=length(X);
M=length(Y);
if(N~=M)
disp('数据点坐标不匹配！');
return;
end
c(1:(2*m+1))=0;
b(1:(m+1))=0;
for j=1:(2*m+1)%求出c和b
for k=1:N
c(j)=c(j)+X(k)^(j-1);
if(j<(m+2))
b(j)=b(j)+Y(k)*X(k)^(j-1);
end
end
end
C(1,:)=c(1:(m+1));
for s=2:(m+1)
C(s,:)=c(s:(m+s));
end
A=b'\C;%直接求解法求出拟合系数
11.LZXEC离散试验数据点的线性最小二乘拟合function[a,b]=LZXEC(x,y)
if(length(x)==length(y))
n=length(x);
else
disp('x和y的维数不相等！');
return;
end%维数检查
A=zeros(2,2);
A(2,2)=n;
B=zeros(2,1);
for i=1:n
A(1,1)=A(1,1)+x(i)*x(i);
A(1,2)=A(1,2)+x(i);
B(1,1)=B(1,1)+x(i)*y(i);
B(2,1)=B(2,1)+y(i);
end
A(2,1)=A(1,2);
s=A\B;
a=s(1);
b=s(2);
12.ZJZXEC离散试验数据点的正交多项式最小二乘拟合function a=ZJZXEC(x,y,m)
if(length(x)==length(y))
n=length(x);
else
disp('x和y的维数不相等！');
return;
end%维数检查
syms v;
d=zeros(1,m+1);
q=zeros(1,m+1);
alpha=zeros(1,m+1);
for k=0:m
px(k+1)=power(v,k);
end%x的幂多项式
B2=[1];
d(1)=n;
for l=1:n
q(1)=q(1)+y(l);
alpha(1)=alpha(1)+x(l);
end
q(1)=q(1)/d(1);
alpha(1)=alpha(1)/d(1);
a(1)=q(1);
B1=[-alpha(1)1];
for l=1:n
d(2)=d(2)+(x(l)-alpha(1))^2;
q(2)=q(2)+y(l)*(x(l)-alpha(1));
alpha(2)=alpha(2)+x(l)*(x(l)-alpha(1))^2; end q(2)=q(2)/d(2);
alpha(2)=alpha(2)/d(2);
a(1)=a(1)+q(2)*(-alpha(1));
a(2)=q(2);
beta=d(2)/d(1);
for i=3:(m+1)
B=zeros(1,i);
B(i)=B1(i-1);
B(i-1)=-alpha(i-1)*B1(i-1)+B1(i-2);
for j=2:i-2
B(j)=-alpha(i-1)*B1(j)+B1(j-1)-beta*B2(j); end
B(1)=-alpha(i-1)*B1(1)-beta*B2(1);
BF=B*transpose(px(1:i));
for l=1:n
Qx=subs(BF,'v',x(l));
d(i)=d(i)+(Qx)^2;
q(i)=q(i)+y(l)*Qx;
alpha(i)=alpha(i)+x(l)*(Qx)^2;
end
alpha(i)=alpha(i)/d(i);
q(i)=q(i)/d(i);
beta=d(i)/d(i-1);
for k=1:i-1
a(k)=a(k)+q(i)*B(k);
end
a(i)=q(i)*B(i);
B2=B1;
B1=B;
end
</b)%此循环找出abs(f(x)-p(x))取最大值的x<>。