用matlab实现多次最佳一致的函数逼近

suffer

1. function [sun]=main(n)
   
2. fplot('1/(x+2)',[-1,1],'r');
   
3. x=ones(n+2,1);
   
4. for j=0:n+1
   
5.     x(j+1)=cos(pi*(n+1-j)/(n+1));
   
6. end
   
7. first=ones(n+2,1);
   
8. f=1./(x+2);  %原函数
   
9. last=first;
   
10. for j=2:n+2
    
11.    last(j)=(-1)*last(j-1);
    
12. end
    
13. A=ones(n+2,n+2);
    
14. A(:,1)=first;
    
15. A(:,n+2)=last;
    
16. for j=2:n+1
    
17.     for t=2:j
    
18.         A(:,j)=x.*A(:,j);
    
19.     end
    
20. end
    
21. e=(1e-15)*first; %精度控制条件
    
22. sun=A\f;
    
23. while (1)   
    
24.     at='';
    
25.     for i=1:n
    
26.         for t=1:i
    
27.             if (t==1)
    
28.                 at=strcat('(',num2str(sun(t+1)),')');
    
29.             elseif (t>1)
    
30.                 xt='t';j=2;
    
31.                 while j<t
    
32.                     xt=strcat('t*',xt);j=j+1;
    
33.                 end
    
34.                 at=strcat(num2str(t),'*(',num2str(sun(t+1)),')*',xt,'+',at);
    
35.             end
    
36.         end
    
37.     end
    
38.     %以下得到逼近函数
    
39.     ap1=sun(1:n+1,[1]);
    
40.     for i=1:n+1
    
41.         ap(i)=ap1(n+2-i);
    
42.     end
    
43.     yt=strcat('-1/(t+2)^2=',at);
    
44.     [y]=solve(yt,'t');
    
45.     y=numeric(y);
    
46.     %以下得到一组新的交错点组
    
47.     for i=1:n+1
    
48.         if y(i) < 1 & y(i)>-1
    
49.             for j=2:n+2
    
50.                 if y(i)<x(j)&y(i)>x(j-1)
    
51.                     if (1/(x(j-1)+2)-polyval(ap,x(j-1)))*(1/(y(i)+2)-polyval(ap,y(i)))> 0
    
52.                         x(j-1)=y(i);
    
53.                     elseif (1/(x(j-1)+2)-polyval(ap,x(j-1)))*(1/(y(i)+2)-polyval(ap,y(i)))< 0
    
54.                         x(j)=y(i);
    
55.                     end
    
56.                 end
    
57.             end
    
58.         end   
    
59.     end
    
60.     A=ones(n+2,n+2);
    
61.     A(:,1)=first;
    
62.     A(:,n+2)=last;
    
63.     for j=2:n+1
    
64.         for t=2:j
    
65.            A(:,j)=x.*A(:,j);
    
66.        end
    
67.     end
    
68.     f=1./(x+2);
    
69.     sun1=A\f;
    
70.     if(abs(sun1-sun)<e)
    
71.         break;
    
72.     end
    
73.     sun=sun1;
    
74. end
    
75. hold on;
    
76. funcion=poly2sym(ap);
    
77. ezplot(funcion,[-1,1]);
    
78. num=num2str(n);
    
79. legend('原函数曲线',strcat(num,'次逼近函数曲线'));
    
80. title('最佳逼近比较示意图');
    
81. xlabel('x的取值');
    
82. ylabel('f(x)的取值');
    
83. grid on;
    
84. end

复制代码

本程序作者：邹璇