求pi的近似值可以说是比较经典的问题了。在各种软件环境下,用过包括蒙特卡洛等各种方法求过pi的近似值。今天给大家带来通过数值积分的方法来求pi的近似值,并进行一个简单的误差分析
clear;
clc;
m=1000;
n=2*m;
a=0;
b=1;
h=(b-a)/n;
w_x=1;
k=10;
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));
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];
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
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)');
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)
原理和代码都比较简单,就不做解释了,学过数值积分的同学都懂。更多问题可以咨询我的老师宋士苍。邮箱:songshicang@zzu.edu.cn 宋老师是属于计算数学中有限元分析方向的。当然,假如有工科科研的,做出来偏微分方法解不出来的,也可以找我的老师,他会给你提供一个数值解法。