数学软件计算pai的.pptVIP

* Company Logo 数学实验 ? 的 计 算 上海交大数学系 乐经良 知其然,更知其所 以然. -中国先哲 哪里有数,哪里就 有美. - Proclus 乐经良 数学实验 ? 的 计算 乐经良 通过求 的近似值，了解历史上计算 值的一些方法，包括刘徽割圆术，级数展开，数值积分和蒙特卡洛法等。复习微积分中相关知识，比较它们的差异，了解计算方法对提高计算效率的意义。 数学实验 实验目的 乐经良 实际问题 ?―圆周率, 我们十分熟悉的常数. 你也许能写出 ? = 3.1415926535 用Matlab 可以求出? 到几百位 digits(100) vpa(pi) ans = 3.141592653589793238462643383279502884197 16939937510582097494459230781640628620899 8628034825342117068 但你会计算? 的值吗？你又能用几种方法计算？ 如图 AC ? AD ? DC ? AD ? (OC ? OD ) an 2 ? ? ?1 ? 1 ? ( n )2 ? ? 2 ? 4 ? an2 ) ? ? ? ? ? 乐经良 刘徽割园法 ? 从正六边形开始，逐步求边长与 o B C A D 2 2 2 2 2 2 2 2 a an ?1 ? ( 1 an 2 2 2 4 面积 ? 递推法 设边数为6·2n的正多边形边长为 an 乐经良 于是 ? 的值 ? ? 6 ? 2n ?1 ? Sn ?1 ? 3 ? 2n an (刘徽计算到96边形面积，得到? ? 3.141) 用Matlab计算 m文件 function calpi(n) a(1)=1; for i=1:n-1 a(i+1)=sqrt(2-sqrt(4-a(i)^2)); end S=3*2^(n-1)*a(n) 命令窗口输入 format long g Calpi(5) 乐经良 如何提高精度 ? 提高多边形的边数不能 不能完全达到目的 ? 在Matlab文件中解决 符号运算 function calpi1(n) a(1)=sym(1); for i=1:n-1 a(i+1)=sym(sqrt(2-sqrt(4-a(i)^2))); end S=3*2^(n-1)*a(n); %最后进行数值计算，60为数值计算过程中保留 vpa(S,60) 的有效数字 ? x x 5 arctan x ? x ? ? ? ? ? ( ? 1) n ?1 ? 1 ? ? ? ? ? ( ? 1) 乐经良 ? ? 1 1 3 5 1 2 n ? 1 n ?1 ? 4 ? ? 3 3 5 x 2 n ?1 2 n ? 1 ? 1 ? x 2 ? x 4 ? ? ? (?1) n ?1 x 2 n ? 2 ? ? 1 1 ? x 2 利用幂级数计算 积分导出 取 x=1 乐经良 用Matlab计算 创建m文件 calpi 2.m，内容如下： function calpi2(n) S=0; for i=1:n if mod(i,2)==0 S=S-1/(2*i-1); else S=S+1/(2*i-1); end end S=4*S 乐经良 结果如何 ? calpi2(1000) ans = 3.14059265383979 calpi1(10000) ans = 3.14149265359003 calpi1(20000) ans = 3.14154265358982 精度提高很慢! 1 1 ? ( ) ? ( ) ? ? ? ( ? 1) n ? 1 1 2 n ? 1 2 n ? 1 ? ? ( ) ? ( ) ? ? ? ( ? 1) n ?1 1 2 n ?1 2 n ? 1 3 4 arctan ? ar