题目
求函数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