求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])

原理和程式碼都比較簡單，就不做解釋了，學過數值積分的同學都懂。更多問題可以諮詢我的老師宋士蒼。郵箱：[email protected].edu.cn 宋老師是屬於計算數學中有限元分析方向的。當然，假如有工科科研的，做出來偏微分方法解不出來的，也可以找我的老師，他會給你提供一個數值解法。