数值计算方法(2)
数值计算方法第2章2-1节
（2）计算
f
(
a
2
b)
。
（3）若
f
(
a
2
b
)
0
，计算停止；若
f
(
a
2
b
)
f
(a)
0
，用
若
f
(
a
2
b)
f
(b)
0
，以
a
2
b
代替
a
。
a
2
b
代替
b
；
（4）反复执行第二步与第三步，直到区间长缩小到允许误差范围
之内，此时区间中点即可作为所求的近似解。
18
证明方程 x3 3x2 6x 1 0 在区间（0，1）内有唯一的实根，并
在[－1,－0.25]，[0.5，1.25]，[1.25，2]各区间内至少有一个实根。
10
2.1.3 区间二分法
定理 函数f(x)在[a,b]上单调连续，且f(a)f(b)<0， 则方程f(x)=0在区间[a,b]上有且仅有一个实根x*。
二分法的基本思想 将有根的区间二分为两个小区间，然后判断根在那 个小区间，舍去无根的小区间，而把有根的小区间 再一分为二，再判断根属于哪个更小的区间，如此 反复 ，直到求出满足精度要求的近似根。
5
有根区间
介值定理 若函数 f (x) 在[a, b] 连续，且
f (a) f (b) 0 ，则方程 f ( x) 0 在(a,b) 内至
少有一个实根。将[a, b] 称为 f (x) 的有根区间。
6
2.1.2 逐步搜索法
假设f(x)在区间[a,b]内有一
个实根x*,若 b – a较小，则可 在(a,b)上任取一点x0作为初始 近似根。
数值计算方法习题答案(第二版)(绪论)
数值分析（p11页）4 试证：对任给初值x 0,0)a >的牛顿迭代公式112(),0,1,2,......k ak k x x x k +=+= 恒成立下列关系式：2112(1)(,0,1,2,....(2)1,2,......kk k x k x x k x k +-=-=≥=证明：（1）(21122k k k k k kx a x x x x +-⎫⎛-=+==⎪ ⎝⎭（2） 取初值00>x ，显然有0>k x ，对任意0≥k ，a a x a x x a x x k k k k k ≥+⎪⎪⎭⎫ ⎝⎛-=⎪⎪⎭⎫ ⎝⎛+=+2121216 证明：若k x 有n 位有效数字，则n k x -⨯≤-110218， 而()k k k k k x x x x x 288821821-=-⎪⎪⎭⎫⎝⎛+=-+ nnk k x x 2122110215.22104185.28--+⨯=⨯⨯<-∴>≥ 1k x +∴必有2n 位有效数字。
8 解：此题的相对误差限通常有两种解法. ①根据本章中所给出的定理：（设x 的近似数*x 可表示为m n a a a x 10......021*⨯±=，如果*x 具有l 位有效数字，则其相对误差限为()11**1021--⨯≤-l a x x x ，其中1a 为*x 中第一个非零数）则7.21=x ，有两位有效数字，相对误差限为025.010221111=⨯⨯≤--x x e 71.22=x ，有两位有效数字，相对误差限为025.010221122=⨯⨯≤--x x e 3 2.718x =，有两位有效数字，其相对误差限为：00025.010221333=⨯⨯≤--x e x ②第二种方法直接根据相对误差限的定义式求解 对于7.21=x ，0183.01<-e x∴其相对误差限为00678.07.20183.011≈<-x e x 同理对于71.22=x ，有003063.071.20083.022≈<-x e x 对于718.23=x ，有00012.0718.20003.033≈<-x e x备注：（1）两种方法均可得出相对误差限，但第一种是对于所有具有n 位有效数字的近似数都成立的正确结论，故他对误差限的估计偏大，但计算略简单些；而第二种方法给出较好的误差限估计，但计算稍复杂。
数值计算方法》试题集和答案(1_6)2
《计算方法》期中复习试题一、填空题：1、已知3.1)3(,2.1)2(,0.1)1(===f f f ，则用辛普生（辛卜生）公式计算求得⎰≈31_________)(dx x f ,用三点式求得≈')1(f 。
答案：，2、1)3(,2)2(,1)1(==-=f f f ，则过这三点的二次插值多项式中2x 的系数为 ，拉格朗日插值多项式为 。
答案：-1，)2)(1(21)3)(1(2)3)(2(21)(2--------=x x x x x x x L3、近似值*0.231x =关于真值229.0=x 有( 2 )位有效数字；4、设)(x f 可微,求方程)(x f x =的牛顿迭代格式是( )；答案)(1)(1n n n n n x f x f x x x '---=+5、对1)(3++=x x x f ,差商=]3,2,1,0[f ( 1 ),=]4,3,2,1,0[f ( 0 )；6、计算方法主要研究( 截断 )误差和( 舍入 )误差；7、用二分法求非线性方程 f (x )=0在区间(a ,b )内的根时，二分n 次后的误差限为( 12+-n a b )；8、已知f (1)＝2，f (2)＝3，f (4)＝，则二次Newton 插值多项式中x 2系数为( )； 11、 两点式高斯型求积公式⎰1d )(xx f ≈(⎰++-≈1)]3213()3213([21d )(f f x x f )，代数精度为( 5 )；12、 为了使计算32)1(6)1(41310---+-+=x x x y 的乘除法次数尽量地少，应将该表达式改写为11,))64(3(10-=-++=x t t t t y ，为了减少舍入误差，应将表达式19992001-改写为199920012+ 。
13、 用二分法求方程01)(3=-+=x x x f 在区间[0,1]内的根,进行一步后根的所在区间为 ，1 ,进行两步后根的所在区间为 ， 。
【doc】绘制流网的数值计算方法
绘制流网的数值计算方法水利水运科学研究1994年6月tan一是一或d—dx一0式中口为与轴的夹角,和分别为在和轴向的分量.(2)式为描述流线的流线方程,对二维渗流场都适用.在地下水(假定为不可压缩)的渗流问题中,有渗流连续性方程+:0drdV故流线方程(2)式则为全微分方程(=』z(1)(2)(4)或d妇,)一d+d=Vdy—dx(5)得渗流速度与流函数的关系为y一1:型}J由Darcy定律fa一一到薹}㈣【J及(6)式可得'+雾一(—k)啬+k(耄一雾)(8)式中,,及为Darcy渗透系数;,)为水头函数.可见,当渗流域为各向同性域时(一,一=O),流函数,)满足Laplace微分方程;当域中介质为各向异性体时≠,一≠O),流函数不满足Laplace方程.再由(6)和(7)式可得(k警+b)+(k箬+)一o(9)因此,当渗流域由各向同性体组成时,有+等=o㈣,.av.…或V?V一0(11)表明流网中流线与等水头线互为正交,流网为一正交型曲线网络图.反之,在各向异性庠渗流场中,(10)或(n)式均不等于0.即流网流线与等水头线不正交,流网为一斜交型曲线网第1,z期朱岳明等t绘制流同的数值计算方法络图.三,流网的数值计算在渗流域S内,(6)式的两边分别对y及徽分,整理得流函数的Passion型徽分方程+等+警一=o0.0.…再沿面域S的边界曲线r对(6)式分别乘以和n,相加得(12)式解的Neumaun边界条件+=一(13)式中n和,为边界曲线r的外法线方向余弦i渗流速度和,是水头函数的梯度函数,可用有限元法先计算出水头函数在域中的分布,再求得和.如引入单元内水头函数,y)的有限元法插值式,)一∑Ⅳ.(14)则有一一∑∑k警+b警警+警(15)式中,M和m分别为单元结点的水头值,插值基函数和单元结点数固已有水头函数的有限元离散解,和V只是位置坐标的函数.根据变分原理,由(12)和(13)式构成的流函数的解为下述泛函Ⅱ()的驻值点一盯[吉c丢cc警一一f'一nV,)dP(16)因有V=dPⅡ主(VDdS及一JJV,)dS(16)式.q-简化为)一盯[丢(-f-W¨-+,一(17)把(15)式代/k(17)式并在单元内引入流函数的有限元插值(18)式,);∑Ⅳj(18)则在计算域S内,对(17)式进行有限元法的离散运算及单元传导矩阵和左,右端项的集合组装后,可得标准型有限元法线性支配方程水利水运科学研究1994年6月I-K2{}一{Q)(19)式中的{}为流函数的未知结点值向量.~(18)am(19)式,每个单元对系数矩阵[]和右端项{Q}中相应元素的作用分别为b:盯(警警+警警jacz.:=Ⅱ{警[耋ckiaN,+b警]～警[喜ck警+]cz式中为单元的面域}f,j=1.2'..?,m.(19)式的解就为流函数(z,y)的有限元解.有了水头函数,y)和流函数(z,y)的有限元数值解,再经过有限元解后处理子程序的计算,整理后,即可得绘制流网所需的所有技术数据或直接在计算机上绘制出高精度的流网.事实上,(13)式把渗流场中流函数的强制性第一类边界条件(,)一(,y)转化成了第二类自然边界条件,避开了渗出面上流函数事先不易确定的问题,使流网的数值计算大为简化.为了消除系数矩阵[]的奇异性并使解唯一性,在解(19)式前,可根据问题的性质及其物理意义,对某些结点直接引入已知值,如位于浸润线上结点的值为0.0;位于计算域底部不透水面上结点的值为总渗流量q(或取相对值100).四,算例(一)土提中无压渗流场的流网计算'图2为一均质各向同性土堤中无压Fi.P要誊ed印渗流场流网的计算结果,流线与荨水头"ow.ughhomo~eneoandisotop线互为正交,计算域为ABCDEA,上,下游水头分别为6.0和2.0.在求解水头函数的分布及确定浸润线位置时,采用结点虚流量法[3}在求解流函数(,y)的分布时,令浸润线AD上所有结点的值为0.0.劫一盟赴-rdlD第l,2期朱岳明等绘制谎同的数值计算方法(二)闸坝下成层各向异性地基中流网的计算图3为某一水闸下非均质成层各向异性体地基中流网的计算结果.在计算域中,视闸体混凝土为不透水体,厚 2.0rn,深40.0m的棍凝土防渗板桩也为相对不透水体,约40m深的强透水砂砾覆盖层为均质各向异性透水体,其Dare),渗透张量为吡=(m/d)r2.320,01一lo.o1.16J如/s)覆盖层下为水平向层状岩体,其渗透张量为嘲,一[c…m闸坝上,下游水头分别为20.0m和5.0m,计算中令位于闸底板面上结点的=0.0,位于z一一100.0m,z一150.0m和=一100.0m的边界面上的结点满足=舢(为总渗流量). 图3非均质各向异性闸基域中的斜交型流网F.3Tiltedflowztetofeon|inedseepageflowthrou暮hinhomogeneousandanlsotroplc foundationof8sluice如图3所示,在一一40.0m的地基材料分界面上,因受渗流越流量连续性条件的约束,流线与等势线均发生折射;又因材料均为各向异性体,尽管渗透方向与整体坐标和Y轴同向(==O),流线斜交于等水头线,流网为斜交网络图.这种流网难以用徒手勾画, 也不便由模型试验中得到.124水刺水运科学研究1994年6月(三)非均质土石坝中流网的计算图4为一壤土质宽心墙非均质各向同性体土石坝中渗流流网的计算结果.该坝高30.0m,底宽160.Om,心墙底宽96.Om,坝上,下游水头分别为148.Om和130.Om.整个渗流计算域主要由四种透水性不同的材料组成t一区为壤土质宽心墙体,其渗透系数h;1.16×10一m/s}二区为上,下游坝坡区的粗砂砾填筑区,其岛一L16×10～m/s{三区为深选50.Om的细砂砾覆盖层,其=1.16X10-5m/s;四区为覆盖层下的砂质岩层,其h= 1.16×10一m,s.":,\\汹l—?.'.....'..一1160—————————卜—一l伽———图4壤土质宽心墙土石坝中的流网F.4Flowtofseepagethroughanearthdam该坝虽在各子区域内介质为各向同性体,流线与等水头线正交,但在一70.0m和一120.Om等不同透水体的分界面上流线和等水头线均发生折射,使得流网的形态复杂化.显然,流网也不宜用徒手勾画.五,结语(一)本文介绍的方法适用于各种非均质,各向异性域中复杂流网的绘制工作.与传统的几种方法相比,尤其适用于材料分界面处流线和等水头线均折射及各向异性域中斜交流网的绘制.(二)求解漉函数(z,)时.虽基于已求解得水头函数z,)的基础上,但采用的是饲一种网格.(三)在非均质域中,困在相邻材料区的分界线上流线和等水头线均要发生折射,为了提高流网的计算及绘制精度,对这些分界线区单元,应布置得相对密一些.同理?在水力梯度变化较大的有关阻水及排水区,单元也应布置得小些.∞+占●上第l,2期朱岳明等:绘制浇两的数值箅方法12j参考文献1HarrME.GroundwaterandSeepage.LondonlMCGraw—HillBookCompany.19622AahoJ.Finiteelementseepageflownets.IntJforNumericalaMAnalysis MethodsinGeomechanlcs,1984l8t297—3033ZhuYueming,WangRuyun,XuHongbo.Someadaptivetechnlquesfor solutionoffreesurfaceseepageflowthrougharchdamabutments.In±Procof theIntSymposiumonArchDams,1992 NumericalmethodfordrawingflownetofseepageZhuYuemlng(HohaiUniversity)HeJian(ZheiiangProvinceWaterConservancy,~[anagement)ShaoJingdong(ChenduHydroelectricInvestigation,DesignandResearchInstituteofMOE,M~VR) AbstractAtthefirst.thefundamentaltheoryandcharacteristicoftheflownetsof2-D.DarcyS seepageflowproblemareintroduced.Basedonthevariationalprincipleandfiniteelementmethod,thenumericalmethodandthecorrespondingcalculationexpressionsarepresented indetail.Threeillustrativeexamplesaboutthenumerlzationanddrawingoftheflownets ofseepagesthroughembankmentandsluicefoundationaresuccessfullygiven.Themethod hasfoundamentlysolutedtheproblemthatallofthecomplicatedflownetsingeoengineer—ingmaybenumericallydrawnwithquickandhighprecision.Keywords:seepage,flownet,finiteelementmethod,numericalcalculation,tWO,dimensionaIflOW。
数值方法(第2版)答案
C语言编程习题第二章习题2-25.用二分法编程求6x4 -40x2+9=0 的所有实根。
#include <stdio.h>#include <math.h>#define N 10000double A,B,C;double f(double x){return (A*x*x*x*x+B*x*x+C);}void BM(double a,double b,double eps1,double eps2){int k;double x,xe;double valuea = f(a);double valueb = f(b);if (valuea > 0 && valueb > 0 || valuea <0 && valueb < 0) return;printf("Finding root in the range: [%.3lf, %.3lf]\n", a, b);for(k=1;k<=N;k++) {x=(a+b)/2;xe=(b-a)/2;if(fabs(xe)<eps2 || fabs(f(x))<eps1) {printf("The x value is:%g\n",x);printf("f(x)=%g\n\n",f(x));return;}if(f(a)*f(x)<0) b=x;else a=x;}printf("No convergence!\n");}int main(){double a,b,eps1,eps2,step,start;printf("Please input A,B,C:\n");scanf("%lf %lf %lf",&A,&B,&C);printf("Please input a,b, step, eps1,eps2:\n");scanf("%lf %lf %lf %lf %lf",&a,&b,&step,&eps1,&eps2);for (start=a; (start+step) <= b; start += step) { double left = start;double right = start + step;BM(left, right, eps1, eps2);}return 0;}运行：Please input A,B,C:6 -40 9Please input a,b, step, eps1,eps2:-10 10 1 1e-5 1e-5Finding root in the range: [-3.000, -2.000]The x value is:-2.53643f(x)=-0.00124902Finding root in the range: [-1.000, 0.000]The x value is:-0.482857f(x)=0.00012967Finding root in the range: [0.000, 1.000]The x value is:0.482857f(x)=0.00012967Finding root in the range: [2.000, 3.000]The x value is:2.53643f(x)=-0.00124902有时若把判别语句if(fabs(xe)<eps2 || fabs(f(x))<eps1)改为if(fabs(xe)<eps2 && fabs(f(x))<eps1)会提高精度，对同一题运行结果：Finding root in the range: [-3.000, -2.000]The x value is:-2.53644f(x)=-4.26496e-007Finding root in the range: [-1.000, 0.000]The x value is:-0.482861f(x)=-7.3797e-006Finding root in the range: [0.000, 1.000]The x value is:0.482861f(x)=-7.3797e-006Finding root in the range: [2.000, 3.000]The x value is:2.53644f(x)=-4.26496e-007习题2-35. 请用埃特金方法编程求出x=tgx在4.5（弧度）附近的根。
《数值计算方法》试题集和答案(1_6)2
《计算方法》期中复习试题一、填空题:1、已知/（1）= 1°，/⑵=12 /⑶= 1.3,则用辛普生（辛卜生）公式计算求2、/（1）= 7 /⑵=2, /（3） = 11则过这三点的二次插值多项式中疋的系数为 _____ ,拉格朗日插值多项式为 ________________________ o3、近似值^=0 231关于真值A = 0.229有（2 ）位有效数字;4、设/（X ）可微,求方程x = f^的牛顿迭代格式是（）;5、 对/U ）= X 3 + A + 1J 差商/[0,1,2,3] =（ 1/[0,1,2,3,4] =（）；6、 计算方法主要研究（截断）误差和（舍入）误差；7、用二分法求非线性方程f （x ）=0在区间（s,®内的根时，二分力次后的误差限b_a为（诃 ）;8、已知f （l ）=2, f （2）=3, f （4）=,则二次Newton 插值多项式中#系数为（ ）;f/Wdv打⑴血-+ 八2^）]11、两点式高斯型求积公式J 。
八^（Jo22血2厲 ），代数精度为（5 ）;“34 6 y = 10 H ------ 1 --------- ------------12、 为了使计算x-1 （x-1）- （x-1）-的乘除法次数尽量地少，应将V2ooi - 71999"改写为、/^55T+Vi^ 。
得丄"如 答案:，用三点式求得广⑴〜 ____________答案：心+1 =心答案1 一广（占）该表达式改写为 y = 10+（3 + （4 — 6小）人 t =x-l ,为了减少舍入误差，应将表达式13、用二分法求方程/（X）= Q + x_ 1 = 0在区间［0, 1］内的根，进行一步后根的所在区间为，1 ，进行两步后根的所在区间为，。
14、计算积分匚5低肚，取4位有效数字。
用梯形公式计算求得的近似值为—，用辛卜生公式计算求得的近似值为梯形公式的代数精度为1 ,辛卜生公式的代数精度为3 。
《数值计算方法》试题集及答案(1-6) 2
《计算方法》期中复习试题一、填空题：1、已知3.1)3(,2.1)2(,0.1)1(===f f f ，则用辛普生（辛卜生）公式计算求得⎰≈31_________)(dx x f ,用三点式求得≈')1(f 。
答案：2.367，0.252、1)3(,2)2(,1)1(==-=f f f ，则过这三点的二次插值多项式中2x 的系数为 ，拉格朗日插值多项式为 。
答案：-1，)2)(1(21)3)(1(2)3)(2(21)(2--------=x x x x x x x L3、近似值*0.231x =关于真值229.0=x 有( 2 )位有效数字；4、设)(x f 可微,求方程)(x f x =的牛顿迭代格式是( )；答案)(1)(1n n n n n x f x f x x x '---=+5、对1)(3++=x x x f ,差商=]3,2,1,0[f ( 1 ),=]4,3,2,1,0[f ( 0 )；6、计算方法主要研究( 截断 )误差和( 舍入 )误差；7、用二分法求非线性方程 f (x )=0在区间(a ,b )内的根时，二分n 次后的误差限为( 12+-n a b )；8、已知f (1)＝2，f (2)＝3，f (4)＝5.9，则二次Newton 插值多项式中x 2系数为( 0.15 )； 11、 两点式高斯型求积公式⎰1d )(xx f ≈(⎰++-≈1)]3213()3213([21d )(f f x x f )，代数精度为( 5 )；12、 为了使计算32)1(6)1(41310---+-+=x x x y 的乘除法次数尽量地少，应将该表达式改写为11,))64(3(10-=-++=x t t t t y ，为了减少舍入误差，应将表达式19992001-改写为199920012+ 。
13、 用二分法求方程01)(3=-+=x x x f 在区间[0,1]内的根,进行一步后根的所在区间为 0.5，1 ,进行两步后根的所在区间为 0.5，0.75 。
天津大学《数值计算方法》在线作业二答案
A.按模最大
B.按模最小
C.全部
D.任意一个
?
正确答案：B
8. ()是利用函数的值求自变量的值。
A.三次样条插值
B.反插值
C.分段插值
D.爱尔米特插值
?
正确答案：B
9. A.
B.
C.
D.
?
正确答案：B
10.梯形公式是求解常微分方程的()阶方法。
《数值计算方法》在线作业二
一,单选题
1. A. 1
B. 2
C. 0
D. 3
?
正确答案：A
2.设f(-1)=1,f(0)=3,f(2)=4,则抛物插值多项式中x2的系数为()。
A. -0.5
B. 0.5
C. 2
D. -2
?
正确答案：A
3.求解一阶常微分方程初值问题的梯形公式为()步法。
A.多
B. 2
C. 3
A.错误
B.正确
?
正确答案：A
8.高斯-塞德尔迭代法一定比雅可比迭代法收敛快。()
A.错误
B.正确
?
正确答案：A
9. A.错误
B.正确
?
正确答案：A
10.逐次超松弛迭代法是高斯-赛.正确
?
正确答案：B
A. 2
B. 4
C. 3
D. 5
?
正确答案：A
二,判断题
1.梯形方法是一种隐式的多步法。()
A.错误
B.正确
?
正确答案：A
2.求解微分方程初值问题的向后Euler法是隐式方法。()
A.错误
B.正确
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值计算（一）主讲：张森2011-7-9一、矩阵的数值计算相关MATLAB函数提示：二、插值法1、插值有关的MATLAB函数：2、拉格朗日和牛顿插值法（1） 拉格朗日多项式和基函数的MATLAB 程序 求拉格朗日插值多项式和基函数的MATLAB 主程序function [C, L,L1,l]=lagran1(X,Y) m=length(X); L=ones(m,m); for k=1: m V=1;for i=1:m if k~=iV=conv(V,poly(X(i)))/(X(k)-X(i)); endendL1(k,:)=V; l(k,:)=poly2sym (V) endC=Y*L1;L=Y*l例1 给出节点数据03.17)15.2(=-f ，24.7)00.1(=-f ，05.1)01.0(=f ，03.2)02.1(=f ， 06.17)03.2(=f ，05.23)25.3(=f ，作五次拉格朗日插值多项式和基函数，并写出估计其误差的公式.解 在MATLAB 工作窗口输入程序>> X=[-2.15 -1.00 0.01 1.02 2.03 3.25]; Y=[17.03 7.24 1.05 2.03 17.06 23.05]; [C, L ,L1,l]= lagran1(X,Y)运行后输出五次拉格朗日插值多项式L 及其系数向量C ，基函数l 及其系数矩阵L 1如下C =-0.2169 0.0648 2.1076 3.3960 -4.5745 1.0954L =1.0954-4.5745*x+3.3960*x^2+2.1076*x^3+0.0648*x^4-0.2169*x^5 L1 =-0.0056 0.0299 -0.0323 -0.0292 0.0382 -0.0004 0.0331 -0.1377 -0.0503 0.6305 -0.4852 0.0048 -0.0693 0.2184 0.3961 -1.2116 -0.3166 1.0033 0.0687 -0.1469 -0.5398 0.6528 0.9673 -0.0097 -0.0317 0.0358 0.2530 -0.0426 -0.2257 0.0023 0.0049 0.0004 -0.0266 0.0001 0.0220 -0.0002 l =[ -0.0056*x^5+0.0299*x^4-0.0323*x^3-0.0292*x^2+0.0382*x-0.0004] [ 0.0331*x^5-0.1377*x^4-0.0503*x^3+0.6305*x^2-0.4852*x+0.0048] [ -0.0693*x^5+0.2184*x^4+0.3961*x^3-1.2116*x^2-0.3166*x+1.0033] [ 0.0687*x^5-0.1469*x^4-0.5398*x^3+0.6528*x^2+0.9673*x-0.0097] [ -0.0317*x^5+0.0358*x^4+0.2530*x^3-0.0426*x^2-0.2257*x+0.0023] [ 0.0049*x^5+0.0004 *x^4-0.0266*x^3+0.0001*x^2+0.0220*x-0.0002] 估计其误差的公式为)(5x R )25.3)(03.2)(02.1)(01.0()00.1)(15.2(!6)()6(----++=x x x x x x fξ，)3.25,-2.15(∈ξ.拉格朗日插值及其误差估计的MATLAB 程序 拉格朗日插值及其误差估计的MATLAB 主程序function [y,R]=lagranzi(X,Y,x,M) n=length(X); m=length(x); for i=1:mz=x(i);s=0.0; for k=1:np=1.0; q1=1.0; c1=1.0;for j=1:nif j~=kp=p*(z-X(j))/(X(k)-X(j)); endq1=abs(q1*(z-X(j)));c1=c1*j; ends=p*Y(k)+s; endy(i)=s; endR=M*q1/c1;例 2 已知5.030sin = ，1707.045sin = ，0866.060sin = ，用拉格朗日插值及其误差估计的MATLAB 主程序求 40sin 的近似值，并估计其误差.解 在MATLAB 工作窗口输入程序>> x=2*pi/9; M=1; X=[pi/6 ,pi/4, pi/3];Y=[0.5,0.7071,0.8660]; [y,R]=lagranzi(X,Y,x,M)运行后输出插值y 及其误差限R 为y = R =0.6434 8.8610e-004.（2） 牛顿插值法的MATLAB 综合程序求牛顿插值多项式、差商、插值及其误差估计的MATLAB 主程序function [y,R,A,C,L]=newdscg(X,Y,x,M) n=length(X); m=length(x); for t=1:mz=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:nA(i,j)=(A(i,j-1)-A(i-1,j-1))/(X(i)-X(i-j+1));endq1=abs(q1*(z-X(j-1)));c1=c1*j; endC=A(n,n);q1=abs(q1*(z-X(n)));for k=(n-1):-1:1C=conv(C,poly(X(k)));d=length(C);C(d)=C(d)+A(k,k); endy(k)= polyval(C, z); endR=M*q1/c1;L(k,:)=poly2sym(C);例3 将区间 [0,π/2] 分成n 等份)3,2(=n ，用x x f y s i n )(==产生1+n 个节点，求二阶和三阶牛顿插值多项式)(2x P 和)(3x P .用它们分别计算sin (π/7) (取四位有效数字)，并估计其误差.解 首先将名为newdscg.m 的程序保存为M 文件，然后在MA TLAB 工作窗口输入程序>> X1=0:pi/4:pi/2; Y1 =sin(X1); M=1; x=pi/7; X2=0:pi/6:pi/2; Y2 =sin(X2);[y1,R1,A1,C1,P2]=newdscg(X1,Y1,x,M),[y2,R2,A2,C2,P3]=newdscg(X2,Y2,x,M)运行后输出插值y 1=)7/sin()7/(2π≈πP 和y 2=)7/sin()7/(3π≈πP 及其误差限R 1和R 2，二阶和三阶牛顿插值多项式P 2和P 3及其系数向量C 1和C 2，差商的矩阵A 1和A 2如下y1 =0.4548 R1 =0.0282 A1 =0 0 0 0.7071 0.9003 0 1.0000 0.3729 -0.3357 C1 =-0.3357 1.1640 0 P2 =-3024156947890437/9007199254740992*x^2+163820246322191/140737488355328*xy2 =0.4345 R2 =9.3913e-004 A2 =0 0 0 0 0.5000 0.9549 0 0 0.8660 0.6991 -0.2443 0 1.0000 0.2559 -0.4232 -0.1139 C2 =-0.1139 -0.0655 1.0204 0 P3 =-1025666884451963/9007199254740992*x^3-4717668559161127/72057594037927936*x^2+4595602396951547/4503599627370496*x3、高元插值及其MATLAB 程序（1） MESHGRID 命令的功能和调用格式调用格式一 [X,Y] =meshgrid (x,y)例1 已知x =-3:0.2:3；y=x ，计算函数437x z -=e 22y x --的值，并作出函数的图形.解 输入程序>> [X,Y] = meshgrid(-3:.2:3, -3:.2:3); Z =7-3* X.^4 .*exp(-X.^2 - Y.^2), mesh(Z)title(' Z =7-3 X^4 exp(-X^2-Y^2)的图形')运行后输出函数值和图形（略）.例2 作出函数xz +=2e22yx --在区域22≤≤-x ，22≤≤-y 上的图形.解 输入程序>> [X,Y] = meshgrid(-2:.2:2, -2:.2:2);Z = 2+X .* exp(-X.^2 - Y.^2); mesh(Z)title('Z = 2+X exp(-X^2 - Y^2)的图形')运行后输出函数值和图形（略）例3 设节点),,(z y x 中的5:0.5:5- =x ,x y =和函数437x Z -=e 22y x --，作Z 在插值点X ,5:5.0:.93-=Y 5.4:5.0:9.4-=处的双三次插值和二元最近邻插值及其图形.解 （1）双三次插值.输入程序>> [x,y] = meshgrid(-5:0.5:5);z =7-3* x.^4 .* exp(-x.^2 - y.^2); xi=-3.9:0.5:5; yi=-4.9:0.5:4.5; [xi,yi] = meshgrid(xi,yi);zi = interp2(x,y,z,xi,yi, 'cubic'), mesh(xi,yi,zi) hold onplot3(x,y,z,'r.', 'markersize',3*5) hold offxlabel('x'), ylabel('y'), zlabel('z'),title('z =7-3 x^4 exp(-x^2 - y^2) 的双三次插值和数据点的图形')运行后屏幕显示Z 在插值点X ,5:5.0:.93-=Y 5.4:5.0:9.4-=处的双三次插值及其图形(略).（2）二元最近邻插值.输入程序>> [x,y] = meshgrid(-5:0.5:5); z =7-3* x.^4 .*exp(-x.^2 - y.^2);xi=-3.9:0.5:5; yi=-4.9:0.5:4.5; [xi,yi] =meshgrid(xi,yi);zi = interp2(x,y,z,xi,yi, 'nearest'), mesh(xi,yi,zi) hold on,plot3(x,y,z,'r.', 'markersize',3*5), holdoffxlabel('x'), ylabel('y'), zlabel('z')title('z =7-3 x^3 exp(-x^2-y^2) 的二元最近邻插值和数据点的图形')运行后屏幕显示Z 在插值点X ,5:5.0:.93-=Y 5.4:5.0:9.4-=处的二元最近邻插值及其图形(略).（2) 三元插值及其MATLAB 程序例3 设节点),,(z y x 的坐标为y z y x =-=-=,15,3,0,1,12,1,0,4，计算函数x V +=2e 222zy x---在插值点13:25.0:3,3:25.0:3,10:25.0:3-=-=-=i i i z y x 处的三元线性插值，并作其图形.解 输入程序>> x=[-4,0,1,12];y=[-1,0,3,15];z=y; [X,Y,Z]= meshgrid(x,y,z);V= 2+X .* exp(-X.^2 - Y.^2- Z.^2); [xi,yi,zi] =meshgrid(-3:.25:10,-3:.25:3,-3:.25:13);vi = interp3(X,Y,Z,V,xi,yi,zi),slice(xi,yi,zi,vi,[-1 6 9.5],9,[-2 .2 9]), shading flat,lighting flatxlabel('x'), ylabel('y'), zlabel('z'),title('V=2+ x exp(-x^2 - y^2-z^2) 的三元线性插值图形') hold on,colorbar('horiz'), view([-30 45])运行后屏幕显示三元线性插值及其图形（略）.实例：1994年全国数学建模竞赛A 题练习题：（二维插值） 在一丘陵地带测量高程，x 和y 方向每隔100米测一个点，得高程数据如下。