3th_numAna_interpolation 　　

3th_numAna_interpolation

数值分析中的重要内容：插值

26 Sep 2016

Go back 第三次数值分析实验课：插值

一、实验目标

1. 掌握常用的插值方法（Lagrange插值、Newton插值、Hermite插值、三次样条函数插值），应用插值函数求函数在指定点的近似值，并会进行误差估计。

2. 通过实验了解各种插值方法的优缺点。

二、实验问题

1.下列数据表给出了函数\(y=f(x)\)在若干个节点上的函数值\(y_{i}=f(x_{j}),j=(0,1,2,3,4,5)\)。

\(x_{j}\) 0.04 0.55 0.65 0.80 0.95 1.05
\(y_{j}\) 0.41075 0.57815 0.69675 0.90 1.00 1.25382

计算\(f(0.59)\)的值.

1.编写拉格朗日插值函数

以p28的例题为例,利用例题有助于理解,并且可以验证

x<-c(0.32,0.34,0.36)
y<-c(0.314567,0.333487,0.352274)

lagrange<-function(xx,x,y){
  len<-length(x)#插值节点的个数
  l<-vector(mode="numeric",len)#先定义一个向量
  for(i in 1:length(x)){#i用来指示第几个
    others<-setdiff(1:len,i)
    numerator<-1
    denominator<-1
    for (j in others) {
      numerator<-(xx-x[j])*numerator
      denominator<-(x[i]-x[j])*denominator
    }
    l[i]<-y[i]*numerator/denominator
  }
  return(sum(l))
}

lagrange(xx=0.3367,x=x,y=y)
## [1] 0.3303744
#函数正确性通过验证
x<-c(0.04,0.55,0.65,0.80,0.95,1.05)
y<-c(0.41075,0.57815,0.69675,0.90,1.00,1.25382)

利用拉格朗日插值法计算得到的\(f(0.596)\)值为：0.6145127

2.编写牛顿插值函数

以课本p32页为例，利用例题有助于理解,并且可以验证

x<-c(0.40,0.55,0.65,0.80,0.90)
y<-c(0.41075,0.57815,0.69675,0.88811,1.02652)
xx=0.596
#计算均差，用一个矩阵存储
newton<-function(xx,x,y){
  len<-length(x)#插值节点的个数
  mat<-matrix(NA,nrow = len,ncol = len)#存储各阶均差
  mat[,1]<-y
  for (j in 2:len) {#先计算列，j为列数
    for (i in j:len) {
      mat[i,j]<-(mat[i,j-1]-mat[i-1,j-1])/(x[i]-x[i-(j-1)])#分母注意
    }
  }
  # result<-vector()#向量化存储结果，提高效率
  # result[1]<-y[1]
  sum=0;
  for (i in 1:(len)) {
    tmp=1
    if (i==1) {
      sum<-sum+y[1]
    }else{
      for (j in 1:(i-1)) {
        if (j==1) {
          tmp<-tmp*mat[i,i]*(xx-x[j])
        }else{
          tmp<-tmp*(xx-x[j])
        }
        
      }
      sum<-sum+tmp
    }
  }
  return(sum)
}
#注意这里要注意一个问题
newton(xx=0.596,x=x,y=y)#结果正确
## [1] 0.6319175
x<-c(0.04,0.55,0.65,0.80,0.95,1.05)
y<-c(0.41075,0.57815,0.69675,0.90,1.00,1.25382)

利用牛顿插值法计算得到的\(f(0.596)\)值为：0.6215408

  1. 通过对函数\(f(x)=1/(1+25x^{2})\) x属于[-1,1],的高次插值，观察Runge现象。
#runge函数的编写
runge<-function(x){
  1/(1+25*x^2)
}
n<-10
x<-seq(from=-1,to=1,length.out = n+1)
y<-runge(x)

#编写分段样条插值 P44页的题目进行改造
# x<-c(27.7,28,29,30)
# y<-c(4.1,4.3,4.1,3)
# xx=28.5   #测试所用
fenduan<-function(xx,x,y){
  len<-length(x)#插值节点的个数
  u<-0;la<-0;d<-0;h<-0
  #先计算h
  for (i in 1:(len-1)) {
    h[i]<-x[i+1]-x[i]#记住这里下标的1就是书上下标的0，但是这里的h1-3就是书上的h0-2
  }
  h[len]=0#多加一个，为了计算u和la
  #存储各阶均差
  mat<-matrix(NA,nrow = len,ncol = len)
  mat[,1]<-y
  for (j in 2:len) {#先计算列，j为列数
    for (i in j:len) {
      mat[i,j]<-(mat[i,j-1]-mat[i-1,j-1])/(x[i]-x[i-(j-1)])#分母注意
    }
  }
  #计算u la 和 d
  for (i in 2:(len)) {#这里的变量书上是从1开始的，就是说我现在要从2开始
    u[i]=h[i-1]/(h[i-1]+h[i])#u2-4对应课本上的1-3
    la[i]=h[i]/(h[i-1]+h[i])#这里只有la1-(len-1)有意义！
    #注意这里的均差在mat矩阵中的对应关系
    #mat中取值的测试i=2  d[i] mat[i+1,3]
    if(i<len){
      d[i]=6*mat[i+1,3]
    } 
  }
  la[1]<-1#定义la[0]=1
  la[len]<-h[1]/(h[len-1]+h[1])
  u[len]<-1-la[len]
  #利用矩阵形式来进行求解
  A<-matrix(0,nrow = len-1,ncol = len-1)
  for (i in 1:(len-1)) {
      if (i==1) {
        A[i,len-1]=u[2]
        A[i,1]=2
        A[i,2]=la[2]
      }else if(i==(len-1)){
        A[i,1]=la[i+1]
        A[i,(len-2)]<-u[i+1]
        A[i,len-1]<-2
      }else{
        for (j in (i-1):(i+1)) {
          if (j==(i-1)) {
            A[i,j]=u[i+1]
          }else if(j==i){
           A[i,j]=2
          }else{
           A[i,j]=la[i+1]
          }
        }
      }
  }
  #计算所有的M
  d<-d[-c(1)]#第一个元素除去，因为d_{0}不存在，这里多出来了
  d[len-1]<-6*(mat[2,2]-mat[len,2])/(h[len-1]+h[1])
  M=solve(A)%*% d
  xj<-max(which(xx>x))
  xjj<-max(which(xx>x))+1#插值节点的位置
  j<-xj-1
  result<-M[j]*(x[xjj]-xx)^3/(6*h[j+1])+M[j+1]*(xx-x[xj])^3/(6*h[j+1])+(y[xj]-M[j]*h[j+1]^2/6)*(x[xjj]-xx)/h[j+1]+(y[xjj]-M[j+1]*h[j+1]^2/6)*(xx-x[xj])/h[j+1]
  return(result)
}#啊！终于编好了

#为了画图
m=100#取m等于100
x.plot<-seq(from=-1,to=1,length.out = m+1)
y.plot<-runge(x.plot)
y_lagrange<-vector(mode="numeric",length = m+1)
for (i in 1:(m+1)) {
  y_lagrange[i]<-lagrange(xx=x.plot[i],x=x,y=y)
}
#注意fenduan这个函数的rubust有点不行，在这个问题中处理-1~-0.8之间的数据不行
#但是在这个问题中，因为这个runge函数是一个偶函数，所以可以利用偶函数的性质
y_fenduan<-vector(mode="numeric",length = m+1)
for (i in 1:(m+1)) {
  y_fenduan[i]<-fenduan(xx=abs(x.plot[i]),x=x,y=y)
}




#构建一个data.frame，进行画图
df<-data.frame(x=c(x.plot,x.plot,x.plot),y=c(y.plot,y_lagrange,y_fenduan),type=c(rep("f(x)",m+1),rep("P10(x)",m+1),rep("cubic spline interpolation",m+1)))
library(ggplot2)
ggplot(data = df,aes(x=x,y=y,color=type))+geom_point()+geom_line()+geom_hline(aes(yintercept=0),color="yellow")+labs(x="x 的取值",y="利用原函数与插值公式得到的y",title="三个函数式的比较的结果")

从这个两个图可以看出虽然这个插值函数在某些位置的光滑度不是很高，但是我们不可否认的是它对原函数的逼近程度相对于而言还是很好的。而且，当区间分的段数越多三次样条函数就越逼近原函数f(x)。所以我觉得用三次样条插值来解决高次插值出现的龙格现象是很有用的。

3.下面给出误差估计

计算\(P_{10}(x)\)及\(S(x)\)的再点\(x_{i}=-1+pi/12*i,i=1,2...,5\)处的误差，比较两种方法的插值效果

xx<-1:5
xx<- -1+pi/12*xx
#利用拉格朗日插值
y_lagrange<-vector(mode="numeric",length = length(xx))
for (i in 1:length(xx)) {
  y_lagrange[i]<-lagrange(xx=xx[i],x=x,y=y)
}
#利用三次样条插值
y_fenduan<-vector(mode="numeric",length = length(xx))
for (i in 1:length(xx)) {
  y_fenduan[i]<-fenduan(xx=abs(xx[i]),x=x,y=y)
}
#原函数的值
y.plot<-runge(xx)

R中求导估计误差,太费劲了

h<-pi/12
## Higher derivatives:
DD <- function(expr, name, order = 1) {
   if(order < 1) stop("'order' must be >= 1")
   if(order == 1) D(expr, name)
   else DD(D(expr, name), name, order - 1)
}

#计算拉格朗日插值的误差
#得到11阶的导数的情况
eleven1<-DD(expression(1/(1+25*x^2)), "x", 11)
m=50#取m
x.plot<-seq(from=-1,to=1,length.out = m+1)

y<-vector()
for (i in x.plot) {
  x=i
  y<-c(y,eval(eleven1))
}
plot(x.plot,y,type = "l")#可以得到极值为1e15，

#注意：经过我的努力只能进行到这一步了
#求极值
w<-function(xx){
  product<-1
  for (i in 1:length(x)) {
    product<-product*(xx-x[i])
  }
  return(product)
}

m<-1e15;r<-vector()
for (i in 1:length(xx)) {
  r[i]<-m*w(xx[i])/factorial(11)
}

#计算三次样条插值的最大值
tri1<-DD(expression(1/(1+25*x^2)), "x", 4)

m=50#取m
x.plot<-seq(from=-1,to=1,length.out = m+1)

y<-vector()
for (i in x.plot) {
  x=i
  y<-c(y,eval(tri1))
}
plot(x.plot,y,type = "l")

#求极值
tri2<-DD(expression(1/(1+25*x^2)), "x", 5)
library(rootSolve)#载入求根的包
fun <- function (x) {
  eval(as.expression(tri2))
}
All <- uniroot.all(fun, c(-0.3, 0.3))
#得四阶导数的极大值：当x=0的时候导数取得最大值
x=0
m1<-5/384*abs(eval(tri1))*h^4
m2<-1/24*abs(eval(tri1))*h^3
m3<-3/8*abs(eval(tri1))*h^2
error2<-max(m1,m2,m3)#三次样条插值的误差

df<-data.frame(x=c(xx,xx,xx),y=c(y.plot,y_lagrange,y_fenduan),type=c(rep("f(x)",length(xx)),rep("P10(x)",length(xx)),rep("cubic spline interpolation",length(xx))))
#绘图
ggplot(data = df,aes(x=x,y=y,color=type))+geom_point()+geom_line()+geom_hline(aes(yintercept=0),color="yellow")+labs(x="x 的取值",y="利用原函数与插值公式得到的y",title="比较拉格朗日插值和三次样条插值")

最后误差估计的结果

拉格朗日插值估计的结果是-4.35455910^{7}, -3.698696310^{7}, -3.042833710^{7}, -2.38697110^{7}, -1.731108410^{7}

三次样条插值估计的结果是385.5314219

讨论与总结

一、拉格朗日插值与分段低次插值的优缺点： 拉格朗日插值函数在整个区间上有统一的解析表达式，其形式关于节点对称，光滑性也好。但是其收敛性差，稳定性也差。而分段线性插值函数、分段三次埃尔米特插值还有三次样条插值函数虽然光滑性差，但它们都克服了拉格朗日差值函数的缺点，不仅收敛性和稳定性强，而且方法简单实用，计算量也小，所以我觉得它们则更适合被广泛应用于工程应用中，尤其是三次样条插值。因为三次样条插值的光滑度比分段线性插值函数以及分段三次埃尔米特插值都高。

二、不足