当前位置:文档之家› Lagrange插值

Lagrange插值

目录《MATLAB程序设计实践》课程考核 (1)1编程实现“LAGRANGE插值”科学计算算法,并举例应用之 (1)1.1算法说明 (1)1.1.1数学推导 (1)1.1.2Lagrange插值函数 (1)1.2流程图 (3)1.3源代码 (3)1.4运行结果 (4)2编程解决科学计算和工程问题 (5)2.1算法说明 (5)2.2流程图 (6)2.3源代码 (6)2.4运行结果 (7)3求多项式的根并分析误差大小 (8)3.1二分法 (8)3.1.1二分法原理 (8)3.1.2流程图 (9)3.1.3源代码 (9)3.2牛顿迭代法 (10)3.2.1算法说明 (10)3.2.2流程图 (11)3.2.3源代码 (11)3.3以上两种计算方法的运算结果 (12)3.3.1二分法计算结果 (12)3.3.2牛顿迭代法计算结果 (13)《MATLAB 程序设计实践》课程考核1 编程实现“Lagrange 插值”科学计算算法,并举例应用之 1.1 算法说明 1.1.1数学推导由数学理论可知:对12,,,n x x x 个不同的节点,且节点处的取值分别为12,,,n y y y ,则存在插值多项式:()1011n n n P x a a x a x --=+++ , 使得()()1,2,,n i i P x y i n == ,且满足插值条件的次插值多项式唯一。

并基于此定理,推出Lagrange 插值基函数。

考虑一个简单的插值问题:对节点()1,2,i x i n = 中任意一点()1k x k n ≤≤做一n 次多项式()k l x 使它在该点上取值为1,而在其余点()1,2,1,1,,i x i k k n =-+ 上取值为零,即1()0k i i k l x i k=⎧=⎨≠⎩ 表明n 个点()1,2,1,1,i x i k k n =-+ 都是n 次多项式()k l x 的零点,故可设1211()()()()()()k k k k n l x A x x x x x x x x x x -+=----- .其中k A 为待定系数,由条件()1k l x =可得:1111()()()()k k k k k k k n A x x x x x x x x -+=---- .由以上几式联立有:111111()()()()()()()()()k k n k k k k k k k n x x x x x x x x l x x x x x x x x x -+-+----=---- .对应于每一节点()0k x k n ≤≤都能求出一个满足插值条件的n 次插值多项式,这样,由式可以求出1n +个n 次插插多项式12(),(),,()n l x l x l x 。

容易看出,这组多项式仅与节点的取法有关,称它们为在n 个节点上的1n +次基本插值多项式或n 次插值基函数。

1.1.2Lagrange 插值函数利用插值基函数立即可以写出满足插值条件的n 次插值多项式1122()()()n n y l x y l x y l x ++并记为()n L x ;()1122111111111()()()()()()()()()()()()n n n nk k n kk k k k k k k n nn j k k j k jj kL x y l x y l x y l x x x x x x x x x y x x x x x x x x x x y x x -+=-+==≠=++----=-----=-∑∑∏在MATLAB 中不自带Lagrange 函数,需要自行编程实现Lagrange 插值计算。

其功能是利用给出的节点计算出过所有这些节点的插值多项式。

从而达到预测相关工程实践问题中某些物理量的变化趋势。

其算法流程说明如下:(1):,(0,1,2,...,);i i x y i n =输入(2)();L u 计算插值1);u 输入插值点2)0;v =3)0,1,...,k n =对做:o01()/();nk i k i i i kl u x x x =≠=--∏o 2;k k v v l y =+4):,u v 输出。

1.2流程图1.3源代码%lagrange 插值函数%x0,y0为已知的插值点数值%x为所求插值点矩阵%y为返回函数值,即插值函数在所求插值点的函数图1lagrange算法流程图值function 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;%%%%%%%%内循环用于计算y n (x)*l n (x) for i=1:n if i~=jp=p*(z-x0(i))/(x0(j)-x0(i)); end endv=p*y0(j)+v;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end y(k)=v;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% end1.4 运行结果函数调用格式:Lagrange(x0,y0,x)f =。

其中,x0,y0为插值节点,f 是求得的Lagrange 多项式在x 处的函数值。

◆给定ln f x =的数值表1,用上面lagrange 插值计算ln(0.54).>>x0=[0.4:0.1:0.8]; >>y0=log(x0);>>lagrange(x0,y0,0.54)输出结果为:-0.616142*********◆ lagrange 插值中的龙格现象已知函数21()1f x x =+,在区间[-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(x0,y0,'--r')>>hold on>>plot(x0,y1,'-b')图2lagrange 算法中的龙格现象2 编程解决科学计算和工程问题题目:给定由N 个力(1,2,3,,)i F N ⋅⋅⋅组成的平面任意力系,求其合力。

2.1 算法说明◆ 求N 个力的合力,先以矩阵的形式将这N 个力输入到matlab 中储存起来。

因为力是一个矢量,其既有大小也有方向。

因此,我们得分别表示出力的大小以及方向两个参数,并以此为力的特征进行输入。

◆ 任何一个力i F在平面坐标系中都可以对其进行分解,其分解公式为:cos sin ix i iy i F F F F θθ⎧=⋅⎪⎨=⋅⎪⎩ . 其中,θ为i F 与坐标轴X 方向正向的夹角;ix F ,iy F分别为该力在坐标轴,X Y 方向的夹角。

◆ 将所求各个力都分解为两坐标轴方向的后,进行加和,求出两个坐标轴上的合力11n ix x i n iy i F F F F==⎧=⎪⎪⎨⎪=⎪⎩∑∑合合y. ◆ 接着将力合成,力的大小可以表示为:F =合◆ 最后再表示出力与X 轴正方向的夹角即可。

求夹角可以是要注意分类,要根据xF合及F 合y 所在的象限正确的表示出合力与X 轴正向的夹角。

如当0x F > 合,0F >合y 时,可按照下式求解:()arctan yx F F θ=合合.2.2 流程图2.3 源代码 %求合力函数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % f 为个各个力的大小% c 为个各个力与x 轴正向的夹角 % F,C 为所求合力的大小及方向function [F,C]=qiuheli1(f,c)if length(f)~=length(c)disp('输入的力大小的举证与夹角不对应,请重新输入!');a=input('请按顺序输入各个力的大小,并以“[]”括起来,注意中间以空格或逗号隔开:');b=input('请按顺序输入各个力与x 轴正向的图3求合力的流程图夹角,并以“[]”括起来,注意中间以空格或逗号隔开:');qiuheli1(a,b);elsect=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 Fx>0&&Fy>0C=asin(Fy./F)*180/pi;elseif Fx<0&&Fy>0C=asin(abs(Fx)./F)*180/pi+pi/2; elseif Fx<0&&Fy<0C=asin(abs(Fy)./F)*180/pi+pi; elseif Fx>0&&Fy<0;C=asin(abs(Fx)./F)*180/pi+3/2* pi;elseif abs(Fx)==0&&abs(Fy)==0 C=0;elseif Fx>0&&Fy==0C=0;elseif Fx==0&&Fy>0C=90;elseif Fx<0&&Fy==0C=180;elseif Fx==0&&Fy<0C=270;endFCend %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%M文件clc;clear;a=input('请按顺序输入各个力的大小,并以“[]”括起来,注意中间以空格或逗号隔开:');b=input('请按顺序输入各个力与x轴正向的夹角,并以“[]”括起来,注意中间以空格或逗号隔开:');qiuheli1(a,b); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%2.4运行结果调用上面的M文件,然后输入表中的力,可求得合力。

相关主题