prony.m

运用prony法对函数逼近

%prony法函数逼近
clear all
clc
N=35;
M=N;
t=1/N;
for n=1:N
    for m=1:N
        cs(n,m)=f((M-m)*t);
    end
    v(n,1)=-f(M*t);
    M=M+1;
end
a=inv(cs)*v;
a=[1
   a];
z=roots(conj(a'));
s=log(z)/t;

for n=1:N
    for m=1:N
        xs(n,m)=(z(m,1))^(n-1);
    end
    w(n,1)=f((n-1)*t);
end
A=inv(xs)*w;
m=1;
for x=0:0.05:2.5
    r=0;
    for n=1:N
        r=r+A(n,1)*exp(s(n,1)*x);
    end
    y1(m)=r;
    y2(m)=f(x);
    m=m+1;
end
x=0:0.05:2.5;
plot(x,real(y1),'k:');
hold on;
plot(x,y2,'k-');
xlabel('x');
ylabel('f');
legend('prony法','原函数');
title('prony法逼近函数f(x)');
axis tight