声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1686|回复: 0

[共享资源] 用matlab实现多次最佳一致的函数逼近

[复制链接]
发表于 2006-10-10 09:36 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
  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
复制代码


本程序作者:邹璇

评分

1

查看全部评分

回复
分享到:

使用道具 举报

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-12 16:56 , Processed in 0.071766 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表