马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
- function [sun]=main(n)
- fplot('1/(x+2)',[-1,1],'r');
- x=ones(n+2,1);
- for j=0:n+1
- x(j+1)=cos(pi*(n+1-j)/(n+1));
- end
- first=ones(n+2,1);
- f=1./(x+2); %原函数
- last=first;
- for j=2:n+2
- last(j)=(-1)*last(j-1);
- end
- A=ones(n+2,n+2);
- A(:,1)=first;
- A(:,n+2)=last;
- for j=2:n+1
- for t=2:j
- A(:,j)=x.*A(:,j);
- end
- end
- e=(1e-15)*first; %精度控制条件
- sun=A\f;
- while (1)
- at='';
- for i=1:n
- for t=1:i
- if (t==1)
- at=strcat('(',num2str(sun(t+1)),')');
- elseif (t>1)
- xt='t';j=2;
- while j<t
- xt=strcat('t*',xt);j=j+1;
- end
- at=strcat(num2str(t),'*(',num2str(sun(t+1)),')*',xt,'+',at);
- end
- end
- end
- %以下得到逼近函数
- ap1=sun(1:n+1,[1]);
- for i=1:n+1
- ap(i)=ap1(n+2-i);
- end
- yt=strcat('-1/(t+2)^2=',at);
- [y]=solve(yt,'t');
- y=numeric(y);
- %以下得到一组新的交错点组
- for i=1:n+1
- if y(i) < 1 & y(i)>-1
- for j=2:n+2
- if y(i)<x(j)&y(i)>x(j-1)
- if (1/(x(j-1)+2)-polyval(ap,x(j-1)))*(1/(y(i)+2)-polyval(ap,y(i)))> 0
- x(j-1)=y(i);
- elseif (1/(x(j-1)+2)-polyval(ap,x(j-1)))*(1/(y(i)+2)-polyval(ap,y(i)))< 0
- x(j)=y(i);
- end
- end
- end
- end
- end
- A=ones(n+2,n+2);
- A(:,1)=first;
- A(:,n+2)=last;
- for j=2:n+1
- for t=2:j
- A(:,j)=x.*A(:,j);
- end
- end
- f=1./(x+2);
- sun1=A\f;
- if(abs(sun1-sun)<e)
- break;
- end
- sun=sun1;
- end
- hold on;
- funcion=poly2sym(ap);
- ezplot(funcion,[-1,1]);
- num=num2str(n);
- legend('原函数曲线',strcat(num,'次逼近函数曲线'));
- title('最佳逼近比较示意图');
- xlabel('x的取值');
- ylabel('f(x)的取值');
- grid on;
- end
复制代码
本程序作者:邹璇 |