拉格朗日插值多项式积分求圆周率近似Matlab实现

Lagrange 插值多项式积分求圆周率近似
摘要：
公式1：y1=4/(1+x^2) 公式2：y2=4*sqrt(1-x^2) 分别对公式1、公式2求其拉格朗日插值多项式，再对其求0-1上的定积分来求圆周率π的近似值，并在Matlab 中通过画图来比较两个所求得的值与真实值π的偏差。

Lagrange 插值多项式：
)()(l )(L 0i n i n
i x f x x ∑==
其中 )())(())(()
())(())(()(l 11101110i n i i i i i i i n i i x x x x x x x x x x x x x x x x x x x x x -⋯⋯--⋯⋯---⋯⋯--⋯⋯--=+-+-
)(i x f 为函数在i x 处的函数值，)(x L n 为Lagrange 插值多项式。

Matlab 实现：
clc;clear;
a=0;b=1;
n=input('Enter a number n:'); %将0-1分割成n 节点，即n-1段 X=zeros(1,n); %用来放置节点x 的值
P=zeros(1,n); %用来放置节点x 对应的函数值y1 Q=zeros(1,n); %用来放置节点x 对应的函数值y2 x=0;
h=(b-a)/(n-1); %h 为步长
for i=1:n
y1=4/(1+x^2);
y2=4*sqrt(1-x^2);
X(i)=x;
P(i)=y1;
Q(i)=y2;
x=x+h;
end
X;P;Q; %通过循环对X、P、Q进行赋值
syms s;
l=1;z1=0;z2=0;
for j=1:1
for k=2:n
l=l*(s-X(k))/(X(j)-X(k));
end
z1=z1+l*P(j);
z2=z2+l*Q(j);
end
for j=2:n
l=1;
for k=1:j-1
l=l*(s-X(k))/(X(j)-X(k));
end
for k=j+1:n
l=l*(s-X(k))/(X(j)-X(k));
end
z1=z1+l*P(j); %通过循环求的函数y1的Lagrange插值多项式z1 z2=z2+l*Q(j); %通过循环求的函数y2的Lagrange插值多项式z2 end
I1=int(z1,s,0,1); % z1对s在0-1上求定积分
I1=eval(I1) %用小数形式表示I1
I2=int(z2,s,0,1); % z2对s在0-1上求定积分
I2=eval(I2) %用小数形式表示I2
x=3.10:0.0001:3.20;
y0=pi;
y1=I1;
y2=I2;
plot(x,y0,'r') %红线为圆周率π的真实值
hold on
plot(x,y1,'g') %绿线为公式1所求值
hold on
plot(x,y2,'b') %蓝线为公式2所求值
运行结果：
从图中可以看出，当n=6时绿线很接近红线即圆周率π的真实值，而蓝线则偏离较远，当n=11时，绿线基本与红线重叠，而蓝线相对之前来说也减小偏差。

通过比较知道公式1比公式2能更好的求得圆周率π的近似值；随着n的增大，即节点数的增加，Lagrange插值多项式积分所求的值也越来越接近圆周率π的真实值，但是同时也会增大它的运算量，故取一个合适的n使得其结果尽量准确又不至于运算量过大。