Remez算法的MATLAB实现

题目

求函数f(x)=e^x在[-1,1]上的二次多项式逼近。迭代始点为:

(1) x0=-1,x1=-0.5,x2=0.5,x3=1

(2)x0=-1,x1=0,x2=0.5,x3=1

(3)x0=-0.8,x1=-0.3,x2=0.4,x3=0.7

代码实现

利用Remez算法求解。选取初始点,代入方程求解,获得近似多项式和近似偏差。求解近似多项式和函数差值的最小值, 将最小值和近似偏差作差与精度条件比较,若满足条件,则终止迭代,否则利用单一交换法从偏差集中选点替代原点列,继续迭代。

Remez算法介绍:https://wenku.baidu.com/view/784d00717fd5360cba1adbe8.html

function [] = remez(y)
w = 10^(-6);
r4 = 1;
o = 0;
while r4 >= w
    o=o+1;
z = inp(y);
[s,r4] = comp(z);
if r4 < w
    break;
end
y = subs(y,s,z);
end
o
fprintf('%f*x^2+%f*x+%f',z(4),z(3),z(2))
end
 
 
 
function [ x ] = inp( y )
%UNTITLED 此处显示有关此函数的摘要
%   此处显示详细说明
a = zeros(4,4);
b = zeros(4,1);
for i =1:4
    a(i,1) = (-1)^(i-1);
    a(i,2) = 1;
    a(i,3) = y(i);
    a(i,4) = y(i)^(2);
    b(i,1) = exp(y(i));
end
x = pinv(a)*b;
end
 
function [s,r4] = comp(q)
f1 = @(x)exp(x)-q(2)-q(3)*x-q(4)*x^2;
f2 = @(x)q(4)*x^2+q(2)+q(3)*x-exp(x);
[s1,r1] = fminbnd(f1,-1,1);
[s2,r2] = fminbnd(f2,-1,1);
if abs(r1)<=abs(r2)
    r3 = abs(r2);
else
    r3 = abs(r1);
end
s = [s1,s2];
r4 = r3 - abs(q(1));
end
 
function [y] = subs(y,s,q)     %替代,选取新点组,单一交换法
p = length(s);
for i = 1:p
    if s(i)>y(1)||s(i)<y(4)
    for j = 2:4
        if s(i)<y(j)
            c1 = sign(exp(y(j-1))-q(2)-q(3)*y(j-1)-q(4)*(y(j-1))^2);
            c2 = sign(exp(s(i))-q(2)-q(3)*s(i)-q(4)*(s(i))^2);
            if c1 == c2
                y(j-1) = s(i);
            else
                y(j) = s(i);
            end
        end
    end
      else if s(i)> -1 && s(i)<y(1)
            c1 = sign(exp(y(1))-q(2)-q(3)*y(1)-q(4)*y(1)^2);
            c2 = sign(exp(s(i))-q(2)-q(3)*s(i)-q(4)*(s(i))^2);
            if c1 == c2
                y(1) = s(i);
            else
                y(4) = s(i);
            end       
      else if s(i)> y(4) && s(i)<1
            c1 = sign(exp(y(4))-q(2)-q(3)*y(4)-q(4)*y(4)^2);
            c2 = sign(exp(s(i))-q(2)-q(3)*s(i)-q(4)*(s(i))^2);
            if c1 == c2
                y(4) = s(i);
            else
                y(1) = s(i);
            end       
          end 
          end
    end
end
end

计算结果:

1.初始点为-1,-0.5,0.5,1

>> y=[-1,-0.5,0.5,1];
>> remez(y)

0.552894*x^2+1.129678*x+0.989359
r4 =

   7.2588e-07

2.初始点为-1,0,0.5,1

>> y=[-1,0,0.5,1];
>> remez(y)

0.571094*x^2+1.137338*x+0.984055
r4 =

   9.7114e-07

3.初始点为-0.8,-0.3,0.4,0.7

>> y=[-0.8,-0.3,0.4,0.7];
>> remez(y)

0.505991*x^2+1.070909*x+0.999986
r4 =

   4.4691e-07