计算方法课程设计

计算方法课程设计

题目:用Newton插值多项式

处理磁化曲线

学院：理学院

班级：数学 15-1

学生姓名：陈昌林

学生学号： 2014028281

指导教师：石瑞银

2017年 06 月 19 日

课程设计任务书

目录

摘要........................................................................ 错误！未定义书签。

第1章前言............................................................ 错误！未定义书签。

1.1 牛顿差值多项式算法在计算方法课程中地位错误！未定义书签。

1.2牛顿差值多项式算法的发展状况和具体应用领域 (1)

1.3本文主要研究思路与结构安排 (1)

第2章牛顿差值多项式算法基本原理及Matlab程序 (1)

2.1牛顿差值多项式算法的基本原理 (1)

2.2 牛顿差值多项式算法的构造方法 (3)

2.3牛顿差值多项式算法的Matlab程序..... . (4)

第3章利用牛顿差值多项式算法解决实际问题（处理磁化曲线）错误！未定义书签。

3.1 问题提出.................................................. 错误！未定义书签。

3.2 问题的分析与模型建立.......................... 错误！未定义书签。

3.3 数据处理 (7)

3.4 结果分析与误差估计 (9)

结论 (11)

参考文献 (11)

附录 (12)

摘要

在我们生活中，许多实际问题都都可以利用函数来表示某种内在规律的数量关系，但是在很多应用领域，函数)(x f 有时不能直接写出表达式或者过于复杂不易计算，而只能给出函数)(x f 在若干个点上的函数值或导数值，通常也是造一张函数表，当遇到要求表中未列出的变量的函数值时，就必须做数值逼近．例如给定了函数)(x f 在几个特定点上的函数值，为了研究函数的变化规律，就要根据这个表，将其公式化寻求某一函数)(x y 去逼近)(x f ，在给定点上等于函数值，在其他点上约等于函数值，这样)(x y 既能反映)(x f 的函数特性，又便于计算．称)(x y 的插值函数，)(x f 为被插值函数.

我们可以求一个待定函数来近似反映函数的特性，使得待定函数在给定点上等于函数值，在其它点上的函数的值作为函数的近似值，这种方法称为插值法.利用差值函数很容易得到Lagrange 插值多项式，公式结构紧凑，在理论分析中甚为方便.由于公式中的每一项与所有的插值结点有关.因此，如果需要增加一个插值结点，则Lagrange 插值公式中的每一项都需要改变计算量大，为了克服这一缺点，于是产生了Newton 插值法.有的要求更高需要在某些点处的导数也相同于是前两种方法达不到精度于是产生了Hermite 插值法.

本文只讨论牛顿插值法，通过均差的推导引出差值公式，建立Matlab 程序，意在用牛顿插值法处理磁化曲线。

关键词：插值函数，牛顿插值，均差，Matlab 程序

第1章前言

1.1牛顿差值多项式算法在计算方法课程中地位

Newton 插值法是数值逼近中的一个重要部分,它向前继承了Lagrange 插值法,向后引出了Hermite 插值法,可以看作对多项式插值作了一个简单的统一。Newton 插值公式具有形式简单,便于计算等优点。因此,在插值中得到广泛的应用。Newton 插值公式余项更具有一般性，它对于列表函数或)(x f 导数不存在的情形也同样适用。

1.2牛顿差值多项式算法的发展状况和具体应用领域

插值法是一种古老的数学方法，他来自生产实践，早在一千多年前，我国科学家在研究历法上就应用了线性插值与二次插值，但它的基本理论却是在微积分产生之后才逐渐完善的，其应用也逐步增多。特别是在计算机软件中，许三多库函数.等的计算实际上归结于它的逼近函数的计算。在工程的实际问题中，由于航空、造船、精密机械加工等实际问题的需要，插值方法在实践上和理论上显得更为重要，并得到了进一步的发展和广泛的应用。因此我们希望能够得到一个“简单函数”逼近被计算函数的函数值.于是就得到了一种方法。叫做插值逼近或者插值法.

1.3本文主要研究思路与结构安排

文思路是通过均差的定义及性质，导出牛顿插值公式及其误差公式，建立关于牛顿差值的Matlab 程序。在通过具体实际数学问题，用牛顿插值法磁化曲线和铁损曲线。通过本文的研究既加深了对Newton 插值法的理解，同时又了解它的的优缺点。

第2章 牛顿差值多项式算法基本原理及Matlab 程序

2.1牛顿差值多项式算法的基本原理

2.1.1差商的定义：

设函数f

在互异节点 ,,10x x 上的值为

)(0x f , )(1x f ,等，定义 (1)

f

在j i x x ,上的1阶均差为j

i j i j i x x x f x f x x f --=

)

()(],[

(2)f 在k j i x x x ,,上的2阶均差为k

i k j j i k j i x x x x f x x f x x x f --=]

,[],[],,[

(3) 递推地，f 在k x x x ,,,10 上的k 阶均差为：

k

k k k x x x x x f x x x f x x x f --=-02111010]

,,,[],,,[],,,[ 同时规定f 在i

x 上的零阶均差为)

(][]i x f x f =

2.2.2差商的基本性质

性质1 k 阶均差可以表示成1+k

个函数值的线性组合，即

∑

=+-----=k

j k j j j j j j j k x x x x x x x x x f x x x f 011

010)())(()()

(],,,[ 或记为∑=+=k j j k j k x x f x x x f 01

10)(')

(],,,[ω

证明：用数学归纳法。当1=k 时由均差定义有

11100101010)

()()()(],[x x x f x x x f x x x f x f x x f -+-=--=

故上式成立。现假设1-=m k 时已成立，对m k =由均差定义及归纳假设有

m

m m m x x x x x f x x x f x x x f --=-02111010]

,,,[],,,[],,,[

m

m

j m j j j j j j j m

m j m j j j j j j j x x x x x x x x x x x f x x x x x x x x x x x f -----------=∑

∑

=+--=-+-0111

1010111

01))()(()()(1))()(()()(

0121011111

210010201001)())(()(1

)())(())(()

1

1)(

(1

)())(()(x x x x x x x x x f x x x x x x x x x x x x x x x x x f x x x x x x x x x f m m m m m m m

m j m j j j j j j j m j j j m

m -?

---+-?--------+----=--=-+--∑ ∑

=+------=m j m j j j j j j j j x x x x x x x x x x x f 011

10)())(())(()

( 性质2均差对其节点是对称的，即节点按任意顺序排列，其值不变。

性质3如果],,,,[10k x x x x f 是x 的m 次多项式，则其1阶均差],,,,,[110+k k x x x x x f 是x 的1-m 次多项式，且由此递推可得n 阶均差为零阶均差，1+n 阶均差为零。

证明：按均差定义，],,,,[10k x x x x f 的1阶均差为

1

1100110]

,,,[],,,[],,,,[+++--=k k k k x x x x x f x x x f x x x x f

由假设，上式等号右端分子为x 的m 次多项式，且当1+=k x x 时为零，可知分子会有因子)(1+-k x x ,它与分母)(1+-k x x 同时约去，则得等号右端为1-m 次多项式。 j x )(j x f 零阶均差 1阶均差 2阶均差 3阶均差 0x )(0x f ][0x f

1x )(1x f ][1x f ],[10x x f

2x )(2x f ][2x f ],[21x x f ],,[210x x x f

3x )(3x f ][3x f ],[32x x f ],,[321x x x f ],,,[3210x x x x f

2.2.3Newton 差值基本原理

在两点直线公式()()),(,,11++k k k k y x y x 中有：

)()(111k k

k k

k k x x x x y y y x L ---+=++（点斜式）

11111)(++++--+--=k k k k

k k k k y x x x x y x x x x x L （两点式）

我们考虑点斜式，两点为)),(),,((1100y x y x ，则直线方程为：

)()(00

10101x x x x y

y y x N ---+=

那么，在此基础上增加一个节点),(22y x ，则过这三点的插值多项式就是：

)()()(12x c x N x N +=

)(x c 应该是一个二次多项式。根据插值条件

??

?====111120

0102)()()()(y x N x N y x N x N 有0)()(10==x c x c ，所以))(()(10x x x x c x c --= 根据插值条件：

222)(y x N =

可以求出：

)

)(()

())(()()(021221202121122x x x x x N y x x x x x N x N c ---=

---=

重新写)(2x N ：

)()()(12x c x N x N +=

))(()

)(()

()(101202212001010x x x x x x x x x N y x x x x y y y -----+---+

= ))(()(102010x x x x c x x c c --+-+=

其中

)

)(()(120221220

101100x x x x x N y c x x y y c y c ---=

--==

设插值节点为i x ，函数值为i f ，n i ,,2,1,0 =

1,,2,1,0,1-=-=+n i x x h i i i

}max{i h h =

插值条件为n i f x N i i ,,2,1,0,)( == 设插值多项式)(x N 具有如下形式

+--+-+=))(()()(102010x x x x c x x c c x N

)())((110----+n n x x x x x x c

其中n c c c c ,,,,210 为待定系数

2.2 牛顿差值多项式算法的构造方法

假设有n+1个不同的节点及函数在节点上的值（x 0，y 0），……（x n ，y n ），插值多项式有如下形式： ）（）

）（（）（）（）（n 10n 102010n x -x )(x -x x -x x P x x x x x x -??-+??+-++=αααα

（1）

其中系数i α（i=0,1,2……n ）为特定系数，可由插值样条i i n y x P =）（（i=0,1,2……n ）确定。

根据均差的定义，把x 看成[a,b]上的一点,可得 f(x)= f （0x ）+f[10x x ，]（0x -x ）

f[x,0x ]= f[10x x ，]+f[x,10x x ，] （1x -x ） ……

f[x,0x ，…x 1-n ]= f[x,0x ，…x n ]+ f[x,0x ，…x n ](x-x n )

综合以上式子，把后一式代入前一式，可得到：

f(x)= f （0x ）+f[10x x ，]（0x -x ）+ f[210x x x ，，]（0x -x ）（1x -x ）+…+ f[x,0x ，…x n ]（0x -x ）…(x-x n )+ f[x,0x ，…x n ]n ω= N n （x ）+）（x n R 其中

N n （x ）= f （0x ）+f[10x x ，]（0x -x ）+ f[210x x x ，，]（0x -x ）（1x -x ）+…+ f[x,0x ，…x n ]（0x -x ）…(x-x n )（2）

）

（x n R = f(x)- N n （x ）= f[x,0x ，…x n ]n ω（3） n ω=（0x -x ）…(x-x n )

Newton 插值的系数i α（i=0,1,2……n ）可以用差商表示。一般有 f k =α[k 10x x x ??，] （k=0,1,2，……，n ）（4）

把（4）代入（1）得到满足插值条件N ）（）

（i i n x f x =（i=0,1,2，……n ）的n 次Newton 插值多项式

N n （x ）=f （0x ）+f[10x x ，]（1x -x ）+f[210x x x ，，]（1x -x ）（2x -x ）+……+f[n 10x x x ??，]（1x -x ）（2x -x ）…（1-n x -x ）.

其中插值余项为：

）

（）！

（）

（）（）（）（x 1n f x N -x f x R 1n 1

n n +++==ωξ ξ介于k 10x x x ??，之间。

2.3牛顿差值多项式算法的Matlab 程序

2.3.1牛顿插值的程序实现方法：

第一步：计算[][][][]001012012,,,,,,,n f x f x x f x x x f x x x x 、、、、。

第二步：计算牛顿插值多项式中01[,,,]i f x x x ???011()()()i x x x x x x ---???-，

1,2,,i n = ，得到n 个多项式。

第三步：将第二步得到的n 个多项式相加，得到牛顿插值多项式。 第四步：利用所得到的插值多项式，估算x 取其它值时()f x 的值。

第五步：作出所求多项式在插值结点周围的函数图像。

2.3.2牛顿插值法的程序实现

%% newton.m

% 求Newton插值多项式、差商、插值及其误差估计的Matlab主程序[4]

% 输入的量：X是n+1个节点(x_i,y_i)(i=1,2,...,n+1)横坐标向量，Y是纵坐标向量

% x是以向量形式输入的m个插值点，M在[a,b]上满足|f~(n+1)(x)|≤M

% 注：f~(n+1)(x)表示f(x)的n+1阶导数

% 输出的量：向量y是向量x处的插值，误差限R，n次Newton插值多项式L 及其系数向量C

%差商的矩阵A

Function[y,R,A,C,L] = newton(X,Y,x,M)

n = length(X);

m = length(x);

for t = 1:m

z = x(t);

A = zeros(n , n);

A(: , 1) = Y’;

s = 0.0 ; p = 1.0 ; q1 = 1.0 ; c1 = 1.0 ;

for j = 2 : n

for i = j : n

A(i , j) = (A(i , j-1) - A(i-1 , j-1)) / (X(i) - X(i - j + 1));

end

q1 = abs ( q1* ( z - X ( j - 1 ))) ;

c1 = c1 * j ;

end

C = A(n , n) ; q1 = abs ( q1 * ( z - X ( n ))) ;

for k = ( n - 1 ) : -1 : 1

C = conv ( C , poly ( X ( k ))) ;

d = length ( C ) ;

C ( d ) = C ( d ) + A ( k , k ) ;（详细代码见附录一）

第3章利用牛顿差值多项式算法解决实际问题

3.1问题提出

武钢公司在1998年3月向用户发布了新的晶粒取向磁化钢带(片)的磁化曲线和铁损曲线。对于设计人员来讲,查曲线是一件麻烦的事情,不但要拿尺子打准坐标,还要进行估算。尤其是单位铁损曲线,是每次设计都要查的。这样不仅工作量大,而且容易出错。因此,有必要对磁化曲线和铁损曲线进行一些数值处理,让使用者可以根据硅钢片的型号和磁感应强度,方便地找到磁场强度或单位铁损。为此,请采用了数值分析方法中的牛顿插值法对磁化曲线和铁损曲线作了一些处理,以便在实际应用中取得了较好的效果,大大提高工作效率。

3.2 问题的分析与模型建立

首先设法根据表格中已有的函数值来构造一个简单的函数作函数的近似表达式，然后再用来计算新的点上的函数值作为原函数的近似值。通常可以选用多项式函数作为近似函数，因为多项式具有各阶的导数，求值比较方便。用代数多项式作为工具研究插值问题。表格如下

（表一）

3.3 数据处理

根据磁化曲线和铁损曲线的实际情况，采用分段的二次牛顿插值多项式来处理它，在实际应用中完全可以满足工程精度的要求。首先,对于原始曲线,均匀地采集其插值出节点并输入计算机;然后程序采用三点分段插值的办法来构造二次牛顿插值多项式)(2x N ,并根据此计算出所需数据；然后将结果输出。 源程序如下： #include ″stdio.h ″ #include ″iostream.h ″ #include ″math.h ″ void main() {

int n,i ,j ,k,z;

float y[20] ,x[20] ,d[20] ,q[20] ,p,t ,v char chl;

printf(“请输入节点个数n! \ n n ”) ; scanf(″%d,″&n) ;

printf(″请输入各节点值! \ n ”) ; for(i =0;i < =n- 1;i + +) {

scanf(″%f , %f ,″&q[i] ,&d[i]) ; }

for(k=1;k< =20;k+ +) {

printf(″请输入插入点! \ n ″) ; scanf(″%f ,″&t) ; for(i =0;i < =n-1;i + +) {x[i] =q[i] ;y[i] =d[i] ;} for(i =0;i < =n- 1;i + +) / /根据插入点选择合适的节点 {

if(x[i] - t >0) break; }

z=i-1

for(i =1;i < =2;i + +)

{x[i] =x[z] ;y[i] =y[z] ;z=z+1;}

for(i =1;i < =2;i + +)

{

for(j =2;j > =i;j-- )

{

y[j] =(y[j] - y[j - 1])/ (x[j] - x[j - 1]) ; ∥求插商

}

}

printf(″磁感应强度(T)| 铁损(W/ kg) \ n″) ;

V=x[0] ;

for(j =1;j <21;j + +)/ /输出三节点间插值

{

p=y[2] ;

for(i =1;i > =0;i - - )

P=y[i] +(v-x[i])*p;

printf(″%f| ″%f \ n″,v,p) ;

v=v+0.01;/ /步距为0.01

}

cout <<″退出?y/ n\ n″;

cin>>ch1;

if(ch1= =′y′)break;

}

}

此程序用于处理铁损曲线，若需要处理磁化曲线，只需要改动表头汉字，其余步骤均相同。数据结果见表一，（w/kg）差值计算即为最后求解的结果。

3.4 结果分析与误差估计

从中可以发现,插值结果具有良好的收敛性。这是因为实际曲线在一定区间内具有较好的线性特征,所以分段的二次牛顿插值能够达到足够的精度,完全符合实

际生产的需要。

利用牛顿插值法处理磁化曲线和铁损曲线的优点：

(1)用分段的二次牛顿插值法处理磁化曲线和铁损曲线完全可以满足工程精度的

要求。

(2)此法能大大提高工作效率,把设计人员从繁复的查曲线工作中解放出来。

(3)此法具有普遍性。对于其他国家生产的硅钢片,若提供曲线,可同样照此处理。

结论

目前，牛顿插值法已经运用到了工程上的各个领域，并解决了许多实际工程中遇到的问题，如物体加热时间的分析、计算；加药量自动标定；智能气体体积分数测量；自动确定支持度阈值；漏磁探测；电力系统采样；凸轮曲线的修正设计等。有时根据实际情况也会使用局部牛顿插值法，如运用局部牛顿插值法提高多狭缝自准直仪准确度。

在插值问题中，要求插值多项式通过给定的数据点，但实际上所谓给定的数据本身是有误差的，而且即使插值多项式通过了给定的数据点，在这些给定数据点上的误差很小，但在其他点上的误差可能会很大，这是插值问题的缺点。在实际应用中，可以采用与曲线拟合结合等方法来达到更好的效果。

参考文献

[1]李庆扬，王能超，易大义.数值分析[M].北京:清华大学出版社,2009,194-196.

[2]李国霞，牛顿插值法的应用[J].处理磁化曲线，1999，24(3):1-2.

[3]孙绳武,徐萃薇.计算方法引论（第4版）[M].北京:高等教育出版社,2015,85-91.

[4]王建卫，取中水. MATLAB程序设计[M].北京：电子工业出版社，2004，151-152.

附录

附录一：

Function[y,R,A,C,L] = newton(X,Y,x,M)

n = length(X);

m = length(x);

for t = 1:m

z = x(t);

A = zeros(n , n);

A(: , 1) = Y’;

s = 0.0 ; p = 1.0 ; q1 = 1.0 ; c1 = 1.0 ;

for j = 2 : n

for i = j : n

A(i , j) = (A(i , j-1) - A(i-1 , j-1)) / (X(i) - X(i - j + 1));

end

q1 = abs ( q1* ( z - X ( j - 1 ))) ;

c1 = c1 * j ;

end

C = A(n , n) ; q1 = abs ( q1 * ( z - X ( n ))) ;

for k = ( n - 1 ) : -1 : 1

C = conv ( C , poly ( X ( k ))) ;

d = length ( C ) ;

C ( d ) = C ( d ) + A ( k , k ) ;

end

y ( t ) = polyval ( C , z ) ;

R ( t ) = M * q1 / c1 ;

end

L = poly2sym ( C ) ;

课程设计评阅书