函数插值计算(Matlab程序)

原创 2012年03月30日 18:47:22

%Runge函数插值计算
%参考教材:《数值分析》李乃成,梅立泉,科学出版社
%《计算方法教程》第二版 凌永祥,陈明逵
%Runge函数f(x)=(1+25x^2)^(-1),(-1<x<1)
clear;clc;close all;
format short;
fprintf('输入N:(5,10,15,20)\n');
N=input('');
N
i=0:N;
%format short;
fprintf('插值节点xi');xx=-1+2/N*i
fprintf('插值节点xi处的函数值f(xi)');fx=(1+25*xx.^2).^(-1)
yy=fx;
%----------求Newton插值多项式的各阶差商-------
%注意:数组下表为1...N+1,区别于课本中的0...N
for k=1:N
    for i=N+1:-1:k+1
        yy(i)=(yy(i)-yy(i-1))/(xx(i)-xx(i-k));
    end
end
fprintf('Newton插值多项式的各阶差商');yy
syms x y
%--------用秦九韶算法计算多项式的表达式
y=yy(N+1);
for i=N:-1:1
    y=y*(x-xx(i))+yy(i);
end
fprintf('Newton插值多项式的表达式:');simplify(y)
%----------(自然)三次样条插值多项式-----------
M=zeros(1,N+1);%M(i)对应s(x)在节点x(i)处的二阶导数
h=zeros(1,N);%h(i)=x(i)-x(i-1)
A=zeros(N-1,N-1);%三弯矩方程组的系数矩阵(自然三次样条插值)
b=zeros(1,N-1);a=b;c=b;d=b;
h(1)=xx(2)-xx(1);
for i=1:N-1 %i代表各数组中的第i个元素,对应的y的下标为i-1,h为i,d为i+1
    h(i+1)=xx(i+2)-xx(i+1);
    a(i)=h(i)/(h(i)+h(i+1));%三弯矩方程组的第-1条对角线元素
    c(i)=1-a(i);%三弯矩方程组的第-1条对角线元素
    b(i)=2;%三弯矩方程组的主对角线上的元素
    d(i+1)=6/(h(i)+h(i+1))*((fx(i+2)-fx(i+1))/h(i+1)-(fx(i+1)-fx(i))/h(i));
end
d=d(2:N);
M(1)=0;M(N+1)=0;%自然三次样条插值函数的边界条件
%-----------用追赶法解严格三对角占优阵----------
z=zeros(N-1,1);w=z;%z为中间矩阵,w为方程的解
u=zeros(1,N-1);l=u;%u矩阵U的主对角线元素,l矩阵L的第-1条对角线元素
%---------追过程--------
u(1)=b(1);z(1)=d(1);
for i=2:N-1
    l(i)=a(i)/u(i-1);
    u(i)=b(i)-l(i)*c(i-1);
    z(i)=d(i)-l(i)*z(i-1);
end
%---------赶过程--------
w(N-1)=z(N-1)/u(N-1);
for i=N-2:-1:1
    w(i)=(z(i)-c(i)*w(i+1))/u(i);
end
fprintf('三弯矩方程组的解向量:');w'
M(2:N)=w;
for i=1:N
    syms x y;
    fprintf('在第%d个区间(%5.4f<=x<=%5.4f)内,S(x)的表达式为:\n',i,xx(i),xx(i

+1));
    y=(xx(i+1)-x)^3/(6*h(i))*M(i)+(x-xx(i))^3/(6*h(i))*M(i+1)+(yy(i)-h(i)

^2/6*M(i))*(xx(i+1)-x)/h(i)+(yy(i+1)-h(i)^2/6*M(i+1))*(x-xx(i))/h(i);
    y=simplify(y)
end
%-------计算xk(k=1...99)对应的y(xk),N(xk),S(xk)------
K=99;k=1:K;
xk=-1+0.02*k;yk=1./(1+25*xk.^2);
%-------计算N(x)在对应数据点上的的函数------
Nk=zeros(1,K);
for k=1:K
    Nk(k)=yy(N+1);
    for i=N:-1:1
        Nk(k)=Nk(k)*(xk(k)-xx(i))+yy(i);
    end
end
%-------计算S(x)在对应数据点上的的函数------
Sk=zeros(1,K);s=zeros(1,K);
for k=1:K
    for i=1:N-1
        if xk(k)<=xx(i+1)
            s(k)=i;
            break;      %s为xk(k)所在的插值区间号
        else
            s(k)=N;
        end
    end 
    hh=xx(s(k)+1)-xx(s(k));xRight=xx(s(k)+1)-xk(k);xLeft=xk(k)-xx(s(k));
    Sk(k)=(M(s(k))*xRight^3/6+M(s(k)+1)*xLeft^3/6+(fx(s(k))-M(s(k))*hh^2/6)

*xRight+(fx(s(k)+1)-M(s(k)+1)*hh^2/6)*xLeft)/hh;
end
figure(1);plot(xk,yk,'k-',xk,Nk,'k--',xk,Sk,'k:');title('n=20');
legend('原函数f(x)','Newton插值函数N(x)','三次样条插值函数S(x)');
En=abs(Nk-yk);fprintf('Netwon插值误差');EnMax=max(En)
Es=abs(Sk-yk);fprintf('三次样条插值误差');EsMax=max(Es)

 

 

相关文章推荐

函数的插值方法及matlab程序

  • 2009年09月21日 11:49
  • 590KB
  • 下载

牛顿(Newton)插值及其MATLAB程序

拉格朗日插值的优点是格式整齐和规范,有误差估计方式,它的缺点是没
  • seamanj
  • seamanj
  • 2014年07月09日 11:42
  • 10086

数值分析 牛顿插值MATLAB程序

  • 2010年04月29日 10:20
  • 43KB
  • 下载

Newton迭代法求函数极小值点 Matlab程序

牛顿迭代法 是最优化中的

牛顿插值的matlab程序

  • 2009年10月29日 15:19
  • 7.31MB
  • 下载

建模中常用的30各MATLAB程序和函数

内部数学常数 pi  圆周率 exp(1) 自然对数的底数e i 或j               虚数单位 ...

双线性插值matlab程序.doc

  • 2011年05月09日 18:32
  • 31KB
  • 下载

双线性插值放大图像的MATLAB程序

  • 2011年04月24日 19:43
  • 16KB
  • 下载
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:函数插值计算(Matlab程序)
举报原因:
原因补充:

(最多只允许输入30个字)