matlab二次插值法学习笔记

%二次插值法求最小值
clear all
clc

f=@(x)sin(x);
a1=4;
a2=4.5;
a3=5;
Theta_error=0.000000001;
y1=feval(f,a1);
y2=feval(f,a2);
y3=feval(f,a3);
h=10;

c1=(y3-y1)/(a3-a1);
c2=((y2-y1)/(a2-a1)-c1)/(a2-a3);
ap=0.5*(a1+a3-c1/c2);
yp=feval(f,ap);

X=[a1^2 a1 1;
    a2^2 a2 1;
    a3^2 a3 1];
b=inv(X)*[y1;y2;y3];
j=1;
for i=4:0.05:5
    z(j)=b(1)*i^2+b(2)*i+b(3);
    w(j)=i;
    q(j)=sin(i);
    j=j+1;

end
plot(w,q,'m:')
   hold on
   plot(w,z,'b--')
   hold on

stepNum=0;
k=stepNum;
Da1=[];               %为迭代中的各个变化创建存储空间
Da2=[];
Da3=[];
Dy1=[];
Dy2=[];
Dy3=[];
Dc1=[];
Dc2=[];
Dap=[];
Dyp=[];

while abs((y2-yp)/y2) >= Theta_error
    stepNum=stepNum+1;
    k=stepNum;
    if (ap-a2)*h>0
        if y2>=yp
            a1=a2;  Da1(k)=a1;               %每一次迭代后把变化后的值赋给D
            y1=y2;  Dy1(k)=y1;
            a2=ap;  Da2(k)=a2;
            y2=yp;  Dy2(k)=y2;
        else
            a3=ap;  Da3(k)=a3;
            y3=yp;  Dy3(k)=y3;
        end
    else
        if y2>=yp
            a3=a2;  Da3(k)=a3;
            y3=y2;  Dy3(k)=y3;
            a2=ap;  Da2(k)=a2;
            y2=yp;  Dy2(k)=y2;
        else
            a1=ap;  Da1(k)=a1;
            y1=yp;  Dy1(k)=y1;
        end
    end
    c1=(y3-y1)/(a3-a1);  Dc1(k)=c1;
    c2=((y2-y1)/(a2-a1)-c1)/(a2-a3);  Dc2(k)=c2;
    ap=0.5*(a1+a3-c1/c2);      Dap(k)=ap;
    yp=feval(f,ap);            Dyp(k)=yp;
end

if y2<yp
    aa=a2;
    yy=y2;
else
    aa=ap;
    yy=yp;
end

x=4:0.0001:5;
y=sin(x);
plot(x,y,'k')
hold on
plot(aa,yy,'r*')