1、1实验二、函数的插值法一、插值算法的基本思想：函数的变化规律往往是通过一组实验数据给出，为了研究此变化规律往往需要求出不在表上的一些函数值。因此我们希望根据给定的函数表做一个既能反函数 的特性，又便于计算的简单函数 ，用 近似 。通常选)(xf )(xp)()(xf一类较简单的函数作为 ，并使 这样确定的)(xp成 立 。对 nifi ,10)(就是我们希望得到的插值函数，主要有Lagrange 插值；Newton 插值；)(xpHermit 插值；分段线性插值；三次样条插值。二、实验要求：1、利用 Lagrange 插值公式编写插值多项式程序。2、给出分段二次和三次多项式的一般表达式，并根据
2、节点选取规则对上述两种插值法进行比较分析。3、对此插值问题若用牛顿插值多项式结果如何？牛顿插值多项式公式如下：三、目的和意义：1、学会常用的插值计算方法以求函数的近似表达式， 以解决其它科学实验的计算问题；2、通过比较加深理解插值多项式与分段插值多项式的优劣问题；3、熟悉各种插值方法的程序设计；4、通过画出的各种插值曲线比较起插值的光滑度问题。四、算法步骤：0()()nkLxly0()1,()(0,12,)nji ijji jillxn10010010(),()(),knn njkjiinkiijjNxffxxff x2六、Newton 插值1）输入各个节点和节点的函数值 ；ymx,10,10
3、 和2）构造均差表，定义 ， ， ；)(xf,nf3）定义来得 ， ， 代入 Newton 公式，计算 N 次插值多项式f0,nxf,0 流程图：开始 ),10(;,niyxi读 入k),10(;)()1( nIxFIfIf 1);(kFky否是 21x；输 出)();(1nFn)0,12,(*211nIIyx输 21;nk3结束七、MATLAB 程序源代 码Newton 插值程序源代码function f=newton(x0,y0) %x0 为入的节点值，y0 相应节点的函数值n=length(x0)syms xfor i=1:nf(i,1)=y0(i) endhx=f(1,1)xx=(1.
4、0)for k=2:nfor i=k:nf(i,k)=(f(i,k-1)-f(i-1,k-1)/(x0(i)-x0(i-k+1) %构造差商表endxx=xx*(x-x0(k-1)hx=hx+f(k,k)*xxendf=expand(hx) %多项式展开 数值例题1）数值插入：x0=0.40 0.55 0.65 0.80 y0=0.41075 0.57815 0.69675 0.88811运行：newton(x0,y0)2）结果：f=-25542915886587131/11258999068426240000+4572354581681722811/4503599627370496000*x
5、-5029019583899/140737488355328*x2+444355163233905/2251799813685248*x3 数值例题： 多项式 在 取 n=6.按等距节点求分段线性插值函数 L(x):在 matlab 命令2xf6窗口中输入如下程序：a=-6;b=6;n=6;f=x2;fenduan(a,b,n,f)4L = -10*X-24, -6*X-8, -2*X, 2*X, 6*X-8, 10*X-24-6 -4 -2 0 2 4 60510152025303540 数值例题% p60clearx=.25 .3 .39 .45 .53;y=.5 .5477 .6245
6、.6708 .728;% 1% dy=1 .6868;% x0=.4 .47;F,a=scyt(x,y)当 x 属于区间0.25,0.3时这个函数为 f =-6.8*(.30-1.*t)3-4.9*(t-.25)3+.26+.95*t当 x 属于区间0.3,0.39时这个函数为 f =-2.7*(.39-1.*t)3-1.9*(t-.30)3+.30+.85*t当 x 属于区间0.39,0.45时这个函数为 f =-2.9*(.45-1.*t)3-2.2*(t-.39)3+.33+.77*t当 x 属于区间0.45,0.53时这个函数为 f =-1.7*(.53-1.*t)3-1.4*(t-.
7、45)3+.35+.71*ta = 0.6325 0.6855 当 x 属于区间0.25,0.3时这个函数为 f =-6.8*(.30-1.*t)3-4.9*(t-.25)3+.26+.95*t当 x 属于区间0.3,0.39时这个函数为 f =-2.7*(.39-1.*t)3-1.9*(t-.30)3+.30+.85*t当 x 属于区间0.39,0.45时这个函数为 f =-2.9*(.45-1.*t)3-2.2*(t-.39)3+.33+.77*t当 x 属于区间0.45,0.53时这个函数为 f =-1.7*(.53-1.*t)3-1.4*(t-.45)3+.35+.71*ta = 0.
8、6325 0.6855 5clearx=.25 .3 .39 .45 .53;y=.5 .5477 .6245 .6708 .728;% 2% x0=.4 .47;% ddy=0 0;F,a=scyt(x,y)当 x 属于区间0.25,0.3时这个函数为 f =-6.3*(t-.25)3+.26+.97*t当 x 属于区间0.3,0.39时这个函数为 f =-3.5*(.39-1.*t)3-1.6*(t-.30)3+.30+.84*t当 x 属于区间0.39,0.45时这个函数为 f =-2.4*(.45-1.*t)3-2.9*(t-.39)3+.32+.77*t当 x 属于区间0.45,0.
9、53时这个函数为 f =-2.1*(.53-1.*t)3+.36+.70*ta =0.6324 0.68556% p56clearx=27.7 28 29 30;y=4.1 4.3 4.1 3;F,a=scyt(x,y)% 1% dy=3.0 -4;% x0=27.7:0.1:30;当 x 属于区间27.7,28时这个函数为 f =-13.*(28.-1.*t)3+.22*(t-28.)3+19.-.53*t当 x 属于区间28,29时这个函数为 f =.66e-1*(29.-1.*t)3+.14*(t-28.)3+12.-.27*t当 x 属于区间29,30时这个函数为 f =.14*(30
10、.-1.*t)3-1.5*(t-29.)3-12.+.56*ta =4.1000 4.2956 4.3357 4.3000 4.2550 4.2144 4.1787 4.1482 4.1234 4.1047 4.0926 4.0875 4.0898 4.1000 4.1167 4.1318 4.1354 4.1173 4.0678 3.9769 3.8346 3.6310 3.3561 3.0000clearn=7;x=linspace(-pi,pi,n)y=sin(x);F,a=scyt(x,y)% 3% x0=0当 x 属于区间-3.1416,-2.0944时这个函数为 f =.15*(t
11、+3.1)3-3.1-.99*t当 x 属于区间-2.0944,-1.0472时这个函数为 f =.15*(-1.0-1.*t)3+.15*(t+2.1)3-1.2当 x 属于区间-1.0472,0 时这个函数为 f =-.15*t3+.32e-16*(t+1.0)3+.99*t-.39e-16当 x 属于区间0,1.0472时这个函数为 f =.32e-16*(1.0-1.*t)3-.15*t3-.39e-16+.99*t当 x 属于区间1.0472,2.0944时这个函数为 f =-.15*(2.1-1.*t)3-.15*(t-1.0)3+1.27当 x 属于区间2.0944,3.1416
12、时这个函数为 f =-.15*(3.1-1.*t)3+3.1-.99*ta=1.8489e-032或采用：function S,y=bkj(x0,y0,p,x,s1,s2)n=length(x0);w=length(x);for i=1:n-1h(i)=x0(i+1)-x0(i);g(i)=y0(i+1)-y0(i);endfor j=2:n-1d(j)=6*(g(j)/h(j)-(g(j-1)/h(j-1)/(h(j-1)+h(j);u(j)=h(j-1)/(h(j-1)+h(j);r(j)=h(j)/(h(j-1)+h(j);endif p=1r(1)=1;u(n)=1;d(1)=6/h(
13、1)*(g(1)/h(1)-s1);d(n)=6/h(n-1)*(s2-g(n-1)/h(n-1);elseif p=2r(1)=0;u(n)=0;d(1)=2*s1;d(n)=2*s2;elser(n)=h(1)/(h(n-1)+h(1);u(n)=1-r(n);d(n)=6*(g(1)/h(1)-(g(n-1)/h(n-1)/(h(n-1)+h(1);endif p=38for i=1:n-2for j=1:n-2if i=jA(i,j)=2;A(i,j+1)=r(i);A(i+1,j)=u(i+1);endendendA(n-1,n-1)=2;A(1,n-1)=u(1);A(n-1,1)
14、=r(n);for i=1:n-1b(i)=d(i+1);endL=Ab;for i=1:n-1M(i+1)=L(i);endM(1)=M(n);elsefor i=1:n-1for j=1:n-1if i=jA(i,j)=2;A(i,j+1)=r(i);A(i+1,j)=u(i+1);endendendA(n,n)=2;M=Ad;endfor a=1:wif (x(a)=x0(1)R=poly2str(sym2poly(M(j)*(x0(j+1)-s)3/(6*h(j)+M(j+1)*(s-x0(j)3/(6*h(j)+(y0(j)-(M(j)*h(j)2)/6)*(x0(j+1)-s)/h
15、(j)+(y0(j+1)-(M(j+1)*h(j)2)/6)*(s-x0(j)/h(j),x);for t=1:length(R)S(a,t)=R(t);endy(a)=M(j)*(x0(j+1)-x(a)3/(6*h(j)+M(j+1)*(x(a)-x0(j)3/(6*h(j)+(y0(j)-(M(j)*h(j)2)/6)*(x0(j+1)-x(a)/h(j)+(y0(j+1)-(M(j+1)*h(j)2)/6)*(x(a)-x0(j)/h(j);break;endendelse9y(a)=0;disp( x(a) ;disp(this number is not exist it this
16、 erea)endend(注释: 第一种边界条件) 0.31429.87)(ixf 0.4)(372SX=27.75 28.85 29.95Y= 4.2222 4.0877 3.1888得到的 S(x)= 13.293 x3 - 1116.4136 x2 + 31253.5671 x - 291635.1316 ： 0.072277 x3 - 5.8733 x2 + 158.4366 x - 1413.9139 )28,7.(x： -1.6574 x3 + 144.6109 x2 - 4205.604 x + 40771.8119 9CPU 时间: 0.0250s结果分析：本程序主要特点是通用性比较长,三种边界条件都可运用,而且可一次性求出多个点的值.比拉格郎日插值,牛顿插值和分段低次插值 它有光滑性,而且很近似原函数值;10
