Lagrange插值

1、目录MATLAB程序设计实践课程考核11编程实现“LAGRANGE插值”科学计算算法，并举例应用之11.1算法说明11.1.1数学推导11.1.2Lagrange插值函数21.2流程图31.3源代码31.4运行结果42编程解决科学计算和工程问题52.1算法说明52.2流程图62.3源代码62.4运行结果73求多项式的根并分析误差大小83.1二分法83.1.1二分法原理83.1.2流程图93.1.3源代码93.2牛顿迭代法103.2.1算法说明103.2.2流程图113.2.3源代码113.3以上两种计算方法的运算结果123.3.1二分法计算结果123.3.2牛顿迭代法计算结果132MATLAB

2、程序设计实践课程考核1 编程实现“Lagrange插值”科学计算算法，并举例应用之1.1 算法说明1.1.1 数学推导由数学理论可知：对个不同的节点，且节点处的取值分别为，则存在插值多项式：，使得，且满足插值条件的次插值多项式唯一。并基于此定理，推出Lagrange插值基函数。考虑一个简单的插值问题：对节点中任意一点做一次多项式使它在该点上取值为1,而在其余点上取值为零， 即表明个点都是次多项式的零点，故可设.其中为待定系数，由条件可得：.由以上几式联立有：.对应于每一节点都能求出一个满足插值条件的次插值多项式，这样，由式可以求出个次插插多项式。容易看出，这组多项式仅与节点的取法有关，称它们为

3、在个节点上的次基本插值多项式或次插值基函数。1.1.2 Lagrange插值函数利用插值基函数立即可以写出满足插值条件的次插值多项式 并记为；在MATLAB中不自带Lagrange函数，需要自行编程实现Lagrange插值计算。其功能是利用给出的节点计算出过所有这些节点的插值多项式。从而达到预测相关工程实践问题中某些物理量的变化趋势。其算法流程说明如下：开始Yes 结束YesNoNo1.2 流程图图 1 lagrange算法流程图1.3 源代码15% lagrange 插值函数% x0,y0为已知的插值点数值% x为所求插值点矩阵% y为返回函数值，即插值函数在所求插值点的函数值functio

4、n y=lagrange(x0,y0,x)n=length(x0);m=length(x);%最外层循环用于输出结果for k=1:m z=x(k); v=0.0; %外循环用于计算求和 for j=1:n p=1.0; %内循环用于计算yn(x)*ln(x) for i=1:n if i=jp=p*(z-x0(i)/(x0(j)-x0(i); end end v=p*y0(j)+v; % end y(k)=v; %end1.4 运行结果函数调用格式：。其中，为插值节点，f是求得的Lagrange多项式在x处的函数值。u 给定的数值表1，用上面lagrange插值计算ln(0.54).表 1

5、函数lnx部分的取值表x0.40.50.60.70.8lnx-0.91629-0.69315-0.51083-0.35777-0.22314输入：x0=0.4:0.1:0.8;y0=log(x0);lagrange(x0,y0,0.54)输出结果为：-0.616142610505419u lagrange插值中的龙格现象已知函数，在区间-5,5之间进行插值，画出插值函数图象，可以比较出：在区间两端插值函数发散比较厉害，精度大大降低。输入：x=-5:1:5;y=1./(1+x.2);x0=-5:0.1:5;y0=lagrange(x,y,x0);y1=1./(1+x0.2);%绘制图象plot(

6、x0,y0,-r)hold onplot(x0,y1,-b)图 2 lagrange算法中的龙格现象2 编程解决科学计算和工程问题 题目：给定由个力组成的平面任意力系，求其合力。2.1 算法说明u 求N个力的合力，先以矩阵的形式将这N个力输入到matlab中储存起来。因为力是一个矢量，其既有大小也有方向。因此，我们得分别表示出力的大小以及方向两个参数，并以此为力的特征进行输入。u 任何一个力在平面坐标系中都可以对其进行分解，其分解公式为：.其中，为与坐标轴方向正向的夹角；，分别为该力在坐标轴方向的夹角。u 将所求各个力都分解为两坐标轴方向的后，进行加和，求出两个坐标轴上的合力.u 接着将力合成

7、，力的大小可以表示为：.u 最后再表示出力与轴正方向的夹角即可。求夹角可以是要注意分类，要根据及所在的象限正确的表示出合力与轴正向的夹角。如当，时，可按照下式求解：.开始将个平面力进行分解及将坐标轴上的力进行合成及根据合力所在的象限求出合力与轴方向 的夹角输出结果结束2.2 流程图图 3 求合力的流程图2.3 源代码%求合力函数% f为个各个力的大小% c为个各个力与x轴正向的夹角% F,C为所求合力的大小及方向function F,C=qiuheli1(f,c)if length(f)=length(c)disp(输入的力大小的举证与夹角不对应，请重新输入！);a=input(请按顺序输入各

8、个力的大小，并以“”括起来，注意中间以空格或逗号隔开：); b=input(请按顺序输入各个力与x轴正向的夹角，并以“”括起来，注意中间以空格或逗号隔开：); qiuheli1(a,b);else ct=c*pi/180; c1=sym(ct); fx=f.*cos(c1); fy=f.*sin(c1); FX=sum(fx,2); FY=sum(fy,2); F0=(FX.2+FY.2).0.5; format long; F=double(F0); Fx=double(FX); Fy=double(FY); if Fx0&Fy0 C=asin(Fy./F)*180/pi; elseif F

9、x0C=asin(abs(Fx)./F)*180/pi+pi/2; elseif Fx0&Fy0&Fy0&Fy=0 C=0; elseif Fx=0&Fy0 C=90; elseif Fx0&Fy=0 C=180; elseif Fx=0&Fy0 C=270; end F Cend%M文件clc;clear;a=input(请按顺序输入各个力的大小，并以“”括起来，注意中间以空格或逗号隔开：);b=input(请按顺序输入各个力与x轴正向的夹角，并以“”括起来，注意中间以空格或逗号隔开：);qiuheli1(a,b);%2.4 运行结果调用上面的M文件，然后输入表中的力，可求得合力。表 2 部

10、分力合成后所得结果序号F1F2F3F4F5F6F合1(30,0)(40,90)(50,53.13)2(10,45)(10,135)(14.14,90)3(10,0)(20,0)(30,0)(60,0)4(10,0)(10,120)(10,240)(0,0)5(10,0)(10,120)(10,240)(10,360)(10,480)(10,600)(0,0)6(10,0)(10,90)(10,180)(10,270)(0,0)7(10,0)(20,120)(30,240)(40,360)(50,480)(60,600)(34.64,33.14)8(20,0)(30,120)(40,240)(5

11、0,360)(60,480)(70,600)(34.64,33.14)9(10,0)(10,120)(10,240)(10,360)(10,0)10(11,42)(12,145)(13,4)(14,415)(15,945)(16,645)(12.89,2.48)其中（a，b）表示力的两个参数，a为大小，b为方向。3 求多项式的根并分析误差大小3.1 二分法3.1.1 二分法原理二分法是方程求解最直观、最简单的方法。二分法以连续函数的介值定理为基础的。由介值定理知道，若函数区间上连续，且，即和负号相反，则在内一定有实根。二分法的基本思想是：用对分区间的方法根据分点处函数的符号逐步将有限区间缩小，

12、使在足够小的区间内，方程有且仅有一根。下面简述其基本步骤。首先记。用中点将区间等分成2个小区间：和。然后分析可能存在的三种情况：u 如果，则是零点，也就是方程的根。u 如果，则区间内存在零点。u 如果，则区间内存在零点。对有根的新区间施行同样的操作，于是得到一系列有空的区间：其中每1个区间的长度都是前一区间长度的一半，最后1个区间的长度为：如果取最后1个区间的中点：作为根的近似值，则有误差估计式：对于所给精度，若取使得则有3.1.2 流程图图 4 二分法流程图开始选定初始区间取区间的中点中点函数值是否为0MN结束否否是是注意：1. 初始区间是一个两端函数值符号相反的区间.2. “M”的意思是取

13、新区间，其中一个端点是原区间端点，另一个端点是原区间的中点。找新区间的法则是：两个端点的函数值异号。3. “N”的意思是方程的解满足要求的精确度.3.1.3 源代码用二分法求解方程在有根区间内的一个根，其中在只有一个根的情形。% % fun是所要求解的函数% a和b分别为有根区间的左右限% eps是允许的误差界% y为所求近似解% yc为函数f在c上的函数值% err是c的误差估计function y,err,yc=erfen(fun,a,b,eps) if nargin4 eps=1e-4; end if feval(fun,a)*feval(fun,b)eps if feval (fun,

14、a)*feval(fun,c)0 b=c;c=(a+b)/2; elseif feval(fun,c)*feval(fun,b)0 a=c;c=(a+b)/2; else y=c;eps=10000; end n=n+1; d=abs(a-b); end format long; err=d yc=feval(fun,c) y=c elseif feval (fun,a)=0 y=a elseif feval(fun,b)=0 y=b; else disp(这个区间内可能不是有根区间，请重新选定区间计算); end% 函数定义function y=f0(x) y=x2-3*x+2;end%fu

15、nction y=f0(x) y=x2-3*x+2;end%3.2 牛顿迭代法3.2.1 算法说明构代迭代函数的一条重要途径是用近似方程来代替原方程去求根。因此，如果能将非线性方程用线性方程去代替，那么，求近似根问题就很容易解决，而且十分方便。牛顿(Newton)法就是一种将非线性方程线化的一种方法。设是方程的一个近似根，把如果在处作一阶Taylor展开，即：于是我们得到如下近似方程：设，则方程的解为：取作为原方程的新近似根，即令：，上式称为牛顿迭代格式。用牛顿迭代格式求方程的根的方法就称为牛顿迭代法，简称牛顿法。牛顿法具有明显的几何意义。方程：是曲线上点处的切线方程。迭代格式就是用切线式的零

16、点来代替曲线的零点。正因为如此，牛顿法也称为切线法。图 5牛顿迭代法流程图开始输入输出结束NONOYESYES3.2.2 流程图3.2.3 源代码%牛顿法解非线性方程% xr是迭代解% f0是所需求解的方程% df是f0的导数% x0迭代的初始点% eps是要求达到的计算精度function xr=Newtonm(x0,eps)format long;f01=feval(f0,x0);df1=feval(df,x0);if df1=0; error(df(x)/dx=0 at x0);endif nargineps; %注意abs(d)=abs(x1-x0) x1=x0-d; x0=x1; f

17、01=feval(f0,x0); df1=feval(df,x0); if df1=0; error(df(x)/dx=0 at x0); end if abs(f01) xr=Newtonm(1)n超过迭代次数，方程可能不存在实根 xr=Newtonm(100)n超过迭代次数，方程可能不存在实根 xr=Newtonm(10000)n超过迭代次数，方程可能不存在实根 xr=Newtonm(-10)n超过迭代次数，方程可能不存在实根上面选用不同的迭代初值，都给出相同的结论，说明该方程可能不存在实数根。 第二题： xr=Newtonm1(10000)err = 6.300276966653822e

18、-005n = 45xr = 1.000000014292257 xr=Newtonm1(0)err = 3.649825436682951e-007n = 71xr = 1.000000364983023 xr=Newtonm1(-10)err = 1.208443495226178e-005n = 30xr = 1.000000000525739 xr=Newtonm1(-10000000)err = 7.005969216603859e-007n = 92xr = 1.000000700598689 xr=Newtonm1(-400)err = 1.658723927627498e-005n = 59xr = 1.000000000990539由以上结果可知1为原方程的解，由于所取的初值不相同，迭代次数也不一样。每次的计算误差也不相同。 第二题： xr=Newtonm2(-0.00001)err = -5.000000000000001e-006n = 0xr = -5.000000000000000e-006 xr=Newtonm2(-0.001)err = -2.499999999992493e-004n = 1xr = -2.49999999999