newton.m

牛顿插值和三次样条插值的matlab实现

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%           计算方法实验            %%%
%%%       ---Newon插值与三次样条插值  %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc                      %清除命令窗的指令
clear all                %清除所有变量
close all                %关闭图像

%计算函数在插值节点的函数值
n=input('输入插值阶次n:');%插值阶次
X=-1:2/n:1;             %产生插值点
disp('f(x)在点xi=-1+2/ni处的函数值f(xi)')
F=1./(1+25*X.^2)        %产生插值节点的函数值

%Newton插值
DQ=zeros(n+1,n+2);     %构造差商表初始赋值
DQ(:,1)=X';            %xi作为差商表的第一列
DQ(:,2)=F';            %f(xi)作为差商表的第二列     
k=0;                   %计算差商表
for j=3:n+2            %逐列计算n阶差商
    k=k+1;
    for i=j-1:n+1      %差商计算 
        DQ(i,j)=(DQ(i,j-1)-DQ(i-1,j-1))/(DQ(i,1)-DQ(i-k,1));
    end
end
DN=diag( DQ(:,2:n+2)); %差商表对角元素
%%%根据差商表对角元素写出Newton插值多项式N(x)%%%
fprintf('N(x)=\n')                      %输出N(x)=
	for i=1:n+1                         %循环输出每一项
        if DN(i)==0                     %若该项的差商值为0则此项不输出
            continue
        end
		fprintf('%2.4f',DN(i))          %输出对应的差商值
		for j=1:i                       %输出每一项的(x-x0)...(x-xn-1)
            fprintf('(x')               %输出(x
			if X(j)<0                   %若节点xi<0,则输出+x(i-1)
				fprintf('+%2.4f',-X(j))
			else
				fprintf('-%2.4f',X(j))  %若节点xi>=0,则xi前加－输出
            end
			fprintf(')')                %输出）
        end
		if i<=n &&DN(i+1)>0             %未输出完且下一项系数>0则输出＋
            fprintf('+')
        end
        if mod(i,5)==0                   %每行输出5项
            fprintf('\n')
        end
    end
X0=-1:2/100:1;                         %将插值区间100等分
Nx=zeros(length(X0),1)';               %存储Newton插值函数值
DetalX=1;                              %Newton多项式在100等分点上的结果
for i=1:n+1
    Nx=DQ(i,i+1)*DetalX+Nx;
    DetalX=DetalX.*(X0-DQ(i,1));
end
plot(X0,Nx,'r')                        %画出Newton插值多项式
xlabel('x轴')
ylabel('y轴')
title('被插函数、Newton插值和三次样条插值')
%三次样条插值
TS=zeros(n+1);               %三弯矩矩阵
h=2/n;                       %区间长度
Mu=1/2; Lambda=1/2;          %等距离hi和hi+1相等，故Mu和Lambda为1/2
%%%假设已知在区间端点的一阶导数值y0,yn %%%
 x=-1:2:1;
 dy=-50./(1+25*x.^2).^2.*x;   %f(x)一阶导数
 dy0=dy(1);dyn=dy(2);         %在端点-1和1处的一阶导数
 TS(1,1)=2;TS(1,2)=1;         %三弯矩阵赋值
 for i=2:n
     TS(i,i-1)=Mu;
     TS(i,i)=2;
     TS(i,i+1)=Lambda;
 end
 TS(n+1,n)=1;TS(n+1,n+1)=2;   
 d=zeros(n+1,1);              %向量d
 d(1)=6/h*((F(2)-F(1))/h-dy0);%三弯矩方程右端向量赋值
 d(n+1)=6/h*(dyn-(F(n+1)-F(n))/h);
 for i=2:n
     d(i)=6*DQ(i+1,4);        %di=6f[x(i-1),x(i),x(i+1)]
 end
     
%%% 利用追赶法求解三弯矩方程 %%% 
TS=[TS,d];                    %三对角增广矩阵
%%%追赶算法%%%
%%%追过程%%%
TS(1,n+2)=TS(1,n+2)/TS(1,1);   %计算y1=d1/l1
for i=1:n
    TS(i,i+1)=TS(i,i+1)/TS(i,i); %ri=ci/li
    TS(i+1,i+1)=TS(i+1,i+1)-TS(i+1,i)*TS(i,i+1);%l(i+1)=b(i+1)-a(i+1)*ri
    TS(i+1,n+2)=(TS(i+1,n+2)-TS(i+1,i)*TS(i,n+2))/TS(i+1,i+1);%y(i+1)=(d(i+1)-a(i+1)*yi)/l(i+1)
end
%%%赶过程%%%
M=zeros(n+1,1);                   %三对角方程组解向量
M(n)=TS(n+1,n+2);                 %Mn=yn
for i=n:-1:1                      %Mi=yi-ri*M(i+1)
    M(i)=TS(i,n+2)-TS(i,i+1)*M(i+1);
end

%%%根据Mi写出三次样条函数插值多项式S(x)%%%
fprintf('\nS(x)=\n')                      %输出S(x)=
	for i=1:n                           %循环输出每个子区间的表达式
        fprintf('%2.4f*[%2.4f(%2.4f-x)^3',1/(6*h),M(i),X(i+1));%输出1/(6h)[(xi-x)^3M(i-1)
        if M(i+1)>=0                    %Mi>0输出＋
            fprintf('+')
        end
        fprintf('%2.4f(x',M(i+1))       %输出Mi(x
        if X(i)>=0                      %x(i-1)>=0输出－x(i-1))^3]
            fprintf('-%2.4f)^3]',X(i))
        else                            %+输出x(i-1))^3]
            fprintf('+%2.4f)^3]',-X(i)) 
        end
        c1=(F(i)-h^2/6*M(i))/h;         %(F(i)-h^2/6*M(i))/h
        if c1>=0                        %(F(i)-h^2/6*M(i))/h>0则输出+
            fprintf('+')
        end
        fprintf('%2.4f(%2.4f-x)',c1,X(i+1));%打印(F(i)-h^2/6*M(i))/h(xi-x)
        c2=(F(i+1)-h^2/6*M(i+1))/h;      %(F(i+1)-h^2/6*M(i+1))/h
        if c2>=0                         %(F(i+1)-h^2/6*M(i+1))/h>0则输出+
            fprintf('+')
        end
        fprintf('%2.4f(x',c2)            %输出(F(i+1)-h^2/6*M(i+1))/h(x
        if X(i)>=0                       %x(i-1)>=0输出-x(i-1)
            fprintf('-%2.4f),',X(i))
        else                             %否则输出+|x(i-1)|
            fprintf('+%2.4f),',-X(i))
        end
       fprintf('     x属于[%2.4f,%2.4f]\n',X(i),X(i+1));%输出x(i-1)并输出子区间范围
    end
    
%%% 三次样条函数的函数值 %%%
Sx=zeros(length(X0),1)';                %存储三次样条函数的函数值
for k=1:length(X0)                      %对101个等分点求函数值
    if k==1                             %对1特殊处理
        i=2;
    else
    i=ceil((k-1)/(100/n))+1;            %将k=100/n*i+1～100/n*(i+1)划分到区间[X(i-1),X(i)]
    end                                 %计算三次样条函数值
     Sx(k)=1/(6*h)*((X(i)-X0(k))^3*M(i-1)+(X0(k)-X(i-1))^3*M(i))+(F(i-1)-h^2/6*M(i-1))*(X(i)-X0(k))/h+(F(i)-h^2/6*M(i))*(X0(k)-X(i-1))/h;
end
hold on;
plot(X0,Sx,'b')                         %画出三次样条插值多项式
 Y=1./(1+25*X0.^2) ;                    %y=f(xk) 
 hold on 
 plot(X0,Y,'m')
if n==5||n==20   
    disp('输出k=1～10,90～99的Xk,Yk,Nk,Sk值：')
    disp('Xk:1～10如下')               %输出X(k),k=1～10         
    Xk1=X0(2:11)
    disp('Xk:90～99如下')              %输出X(k),k=90～99        
     Xk2=X0(91:100)
   disp('Yk:1～10如下')                %输出被插函数函数值Y(k),k=1～10  
    Yk1=Y(2:11)
    disp('Yk:90～99如下')              %输出被插函数函数值Y(k),k=90～99
     Yk2=Y(91:100)
    disp('Nk:1～10如下')               %输出Newton插值函数值N(k),k=1～10 
    Nk1=Nx(2:11)
    disp('Nk:90～99如下')              %输出Newton插值函数值N(k),k=90～99 
     Nk2=Nx(91:100)
    disp('Sk:1～10如下')               %输出三次样条插值函数值S(k),k=1～10
    Sk1=Sx(2:11)
    disp('Sk:90～99如下')              %输出三次样条插值函数值S(k),k=90～99
     Sk2=Sx(91:100)
     
  disp('Newon插值与原函数最大差值为：')  %输出Newton插值与被插函数在指定区间的最大偏差
  EN=max(max(abs(Yk1-Nk1)),max(abs(Yk2-Nk2))) 
  disp('三次样条插值与原函数最大差值为：')%输出三次样条插值与被插函数在指定区间的最大偏差
  ES=max(max(abs(Yk1-Sk1)),max(abs(Yk2-Sk2)))
 end