function-approximation.rar_函数逼近_函数逼近方法

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<b) %此循环找出abs(f(x)-p(x))取最大值的x 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;