matlab数值积分方法求pi的近似值及其比较

求pi的近似值可以说是比较经典的问题了。在各种软件环境下，用过包括蒙特卡洛等各种方法求过pi的近似值。今天给大家带来通过数值积分的方法来求pi的近似值，并进行一个简单的误差分析

clear;
clc;
m=1000;
n=2*m;%设置划分个数n
a=0;%积分上限
b=1;%积分下限
h=(b-a)/n;
w_x=1;%设置高斯积分公式的权值
k=10;%高斯公式使用的正交多项式的数目
%ps:积分函数:编写函数f（x）写入


%% 复化梯形公式求pi
s=0;
for i=1:n-1;
    s=s+f(a+h*i);
end
format long;
Pi_T=4*h*(1/2*f(a)+s+1/2*f(b));
%% 复化辛普森公式求pi
s1=0;s2=0;
for i=1:m-1
    s1=s1+f(a+2*i*h);
    s2=s2+f(a+(2*i-1)*h);
end
s2=s2+f(b-h);
Pi_S=4/3*h*[f(a)+f(b)+2*s1+4*s2];






%% 高斯公式求pi（未复化）
p(1)=sym(a);p(2)=sym(b);
for i=3:k+1
syms x;
delta(i)=int(w_x*x*p(i-1)^2,x,a,b)/int(w_x*p(i-1)^2,x,a,b);
gamma2(i)=int(w_x*p(i-1)^2,x,a,b)/int(w_x*p(i-2)^2,x,a,b);
gamma2(3)=0;
p(i)=(x-delta(i))*p(i-1)-gamma2(i)*p(i-2);
end
p(1:k)=p(2:k+1);
p(k+1)=[];
p%正交多项式集合
%p=expand(p)
%做w(x)=1的k个正交多项式存放到数组p中；
%%
coe=sym2poly(p(k));
xn=sort(roots(coe)');%最高阶正交多项式零点
for i=1:k-1
A(i,1:k-1)=subs(p(i),x,xn);
A=double(A);
end%求权重中的系数矩阵A
B=[int(w_x*p(1)^2,0,1);zeros(k-2,1)];%求权重方程组右端项
w(1:k-1)=inv(A)*B;
w=double(w);%求得权重w（i）
Pi_G=double(4*w*subs(1/(1+x^2),xn)');%高斯公式求得的pi值

%% 复化高斯公式   （ 两点    w（x）=1 ）
s=0;
for i=0:m-1
    s=s+f(a+(2*i+1-sqrt(1/3))*h)+f(a+(2*i+1+sqrt(1/3))*h);
end
Pi_GC=4*h*s;



%% 输出估计值与误差与可视化
  Pi_T
  Pi_S
  Pi_G
  Pi_GC
  error_T=abs(pi-Pi_T)
  error_S=abs( pi-Pi_S)
  error_G=abs(pi-Pi_G)
  error_GC=abs(pi- Pi_GC)
 %plot([error_T error_S  error_G error_GC])

原理和代码都比较简单，就不做解释了，学过数值积分的同学都懂。更多问题可以咨询我的老师宋士苍。邮箱：songshicang@zzu.edu.cn 宋老师是属于计算数学中有限元分析方向的。当然，假如有工科科研的，做出来偏微分方法解不出来的，也可以找我的老师，他会给你提供一个数值解法。