【数值分析实验MATLAB】插值与拟合：拉格朗日插值、多项式插值、分段插值（三次样条插值）、最小二乘法及平方误差

Matlab代码：

format long
x=[11 12 13];
y=[0.190809 0.207912 0.224951];
x1=11.5;
n=length(x)-1;
m=length(x1);
y1=zeros(1,m);
for k=1:m
    z=x1(k);
    y1(k)=0;
    for i=1:n+1
        l=1;
        for j=1:n+1
            if j~=i
                l=l*(z-x(j))/(x(i)-x(j));
            end
        end
        y1(k)=y1(k)+l*y(i);
    end
end
y1=vpa(y1,7);
y1

输出：

y1 =
 
0.1993685

（1）Matlab代码：

n(1)=5;
n(2)=10;
n(3)=20;
l=1;
L=0;

t=-5:0.5:5;
h=1./(1+t.*t);
plot(t,h,'*')
hold on;

syms x;
for k=1:3
   for j=1:n(k)+1
       m(j)=-5+(j-1)*10/n(k);
   end
   for i=1:n(k)+1
       for j=1:n(k)+1
           if(i~=j)
               l=l*(x-m(j))/(m(i)-m(j));
           end
       end
       L=L+f(m(i))*l;
       l=1;
   end
   a=-5:0.5:5;
   %a=-5:0.1:5;
   y=subs(L,x,a);
   plot(a,y);
   hold on;
   L=0;
end

（2）三次样条插值（三弯矩法）

n(1)=5;
n(2)=10;
n(3)=20;

t=-5:0.1:5;
h=1./(1+t.*t);
plot(t,h,'*');
hold on;

syms x
for k=1:3
   for j=1:n(k)+1
   m(j)=-5+(j-1)*10/n(k);
   end
   h=10./n(k);
   for i=1:n(k)+1
    for j=1:n(k)+1
        if(i==j)
            b1(i,j)=4;
        end
        if(abs(i-j)==1)
             b1(i,j)=1;
        end
       if(i~=j&&abs(i-j)~=1)
           b1(i,j)=0;
       end
    end
   end
   for i=1:n(k)+1
       for j=1:n(k)+1
           if(i==j&&j==1)
               b2(i,j)=2;
           end
           if(i==j&&j==n(k)+1)
               b2(i,j)=2;
           end       
       end
   end
   b=b1-b2;
   b=pinv(b);
   d(1)=6.*((f(m(1+1))-f(m(1)))./h-f1(m(1)));
   for i=2:n(k)-1
           d(i)=6.*(f(m(i+1))-2.*f(m(i))+f(m(i-1)))./(h.*h);
   end
   d(n(k)+1)=-6.*((f(m(n(k)+1))-f(m(n(k))))./h-f1(n(k)+1));
   c=b*d';
   syms x
   for i=1:n(k)
      s(i)=f(m(i))+((f(m(i+1))-f(m(i)))./h-(2*c(i)+c(i+1)).*h/6).*(x-m(i))+0.5.*c(i).*(x-m(i))^2+(c(i+1)-c(i)).*(x-m(i))^3./(6.*h);
       a=m(i):h:m(i+1);
      y=subs(s(i),x,a);
      plot(a,y);
      hold on
   end
end

Matlab代码：

format short
%第一题
X=[1 2 3 4 5];
Y=[4 4.5 6 8 8.5];
m=1;
%最小二乘算法
[~,n]=size(X);
S=zeros(1,2*m+1);
T=zeros(m+1,1);
for k=1:2*m+1
    S(k)=sum(X.^(k-1));
end
for k=1:m+1
    T(k)=sum(X.^(k-1).*Y);
end
A=zeros(m+1,m+1);
a=zeros(m+1,1);
for i=1:m+1
    for j=1:m+1
        A(i,j)=S(i+j-1);
    end
end
a=A\T;
a
%平方误差
g=0;
for i=1:n
    y(i)=a(1);
    for j=2:m+1
        y(i)=y(i)+a(j)*X(i)^(j-1);
    end
    g=g+(y(i)-Y(i))^2;
end
fprintf('平方误差为%.4f\n',g);
plot(X,Y,'*');
hold on;
plot(X,y,'r--');
legend('原始数据','拟合曲线');

输出：

a =

    2.4500
    1.2500

平方误差为0.6750

（1）Matlab代码：

%第二题
X=[2 3 4 7 8 10 11 14 16 18 19];
Y=[106.42 108.2 109.5 110 109.93 110.49 110.59 110.6 110.76 111 111.2];
m=2;
%最小二乘算法
[~,n]=size(X);
S=zeros(1,2*m+1);
T=zeros(m+1,1);
for k=1:2*m+1
    S(k)=sum(X.^(k-1));
end
for k=1:m+1
    T(k)=sum(X.^(k-1).*Y);
end
A=zeros(m+1,m+1);
a=zeros(m+1,1);
for i=1:m+1
    for j=1:m+1
        A(i,j)=S(i+j-1);
    end
end
a=A\T;
a
%平方误差
g=0;
for i=1:n
    y(i)=a(1);
    for j=2:m+1
        y(i)=y(i)+a(j)*X(i)^(j-1);
    end
    g=g+(y(i)-Y(i))^2;
end
fprintf('平方误差为%.4f\n',g);
plot(X,Y,'*');
hold on;
plot(X,y,'r--');
legend('原始数据','拟合曲线');

输出：

a =

  106.2927
    0.6264
   -0.0205

平方误差为2.7796

（2）Matlab代码：

function ercf1()
X=[2 3 4 7 8 10 11 14 16 18 19];
Y=[106.42 108.2 109.5 110 109.93 110.49 110.59 110.6 110.76 111 111.2];
X=1./X;
Y=log(Y);
m=1;
%最小二乘算法
[~,n]=size(X);
S=zeros(1,2*m+1);
T=zeros(m+1,1);
for k=1:2*m+1
    S(k)=sum(X.^(k-1));
end
for k=1:m+1
    T(k)=sum(X.^(k-1).*Y);
end
A=zeros(m+1,m+1);
a=zeros(m+1,1);
for i=1:m+1
    for j=1:m+1
        A(i,j)=S(i+j-1);
    end
end
a=A\T;
a(1)=exp(a(1));
X=1./X;
Y=exp(Y);
a
g=0;
for i=1:n
    y(i)=a(1)*(exp(a(2)/X(i)));
    g=g+(y(i)-Y(i))^2;
end
fprintf('平方误差为%.4f\n',g);
plot(X,Y,'*');
hold on;
plot(X,y,'r--');
legend('原始数据','拟合曲线');

输出：

a =

  111.4940
   -0.0903

平方误差为0.4719

ps：运行后的图片未贴出来