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
        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           �s(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           �s(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