matlab第4次课.pptVIP

实验说明 4、用Monte Carlo方法计算圆周率pi。 参考书上102页求pi/4的方法。 * * 综合实验 1.定积分近似计算的Monte Carlo方法 a=0;b=1;m=10000; s=0;H=exp(1);%s设置为落在曲边梯形内的点数 for i=1:m xi=rand();yi=H*rand(); if yiexp(xi) s=s+1; end; %如果随机点落在曲边梯形内,s增加1 end; fprintf(%s%g\n,exp(x)在[0,1]上的积分约等于,H*(b-a)*s/m) 例1 近似计算 2.Matlab近似计算定积分rsums(f,a,b) syms x;rsums(exp(x),0,1) 3.泰勒级数逼近计算器taylortool 4.最大公约数问题 Matlab用于求两个整数的最大公约数的函数是gcd，比如gcd(4,6)将会得到2，数学上有个这样的结论：若d是a,b的最大公约数，则存在整数p,q，使得 pa+qb=d Matlab中可以通过命令[d,p,q]=gcd(a,b)来获得上述的三个数，比如[d,p,q]=gcd(12,32)运行得到 d = 4 p=3 q=-1 可以验证，确实有 3× 12 -1× 32=4 现在考虑这样的一个问题：任取两个正整数，它们互质的可能性有多大？ m=10000;s=0; for i=1:m a=randint(1,2,[1,10^9]); if gcd(a(1),a(2))==1 s=s+1; end end s/m ans = 0.6043 实际上上述所求概率的精确值应该是 证明：设(a,b)=k的概率为pk 在自然数中任取两个数都是偶数的概率显然为1/4,此时a/2与b/2又对应着两个任意的自然数,于是的(a/2 ,b/2)=1概率为p1. a,b均为偶数时, (a ,b)=2的充要条件是(a/2 ,b/2)=1,所以有p2=(1/4)p1,类似地有pk=(1/k2)p1. 由于， 所以 由 所以 syms x; symsum(1/k^2,k,1,inf) 5.勾股数 满足勾股定理 a^2+b^2=c^2的三个正整数a,b,c称为勾股数。 以下程序求出100以内的所有勾股数： l=0; m=[]; for a=1:100 for c=a+1:100 b=sqrt(c^2-a^2); if (b==floor(b))(ba) l=l+1; m(:,l)=[a;b;c]; end end end m 考察c-b=1时的勾股数，将前面程序略加修改可得到下表： 181 145 113 85 61 41 25 13 5 c 180 144 112 84 60 40 24 12 4 b 19 17 15 13 11 9 7 5 3 a 可以证明，此时勾股数的通项可表示为 c-b=2的情况留着自己讨论，现在给出c-b=3时的勾股数，可以得到下表 435 339 255 183 123 75 39 15 c 432 336 252 180 120 72 36 12 b 51 45 39 33 27 21 15 9 a 145 113 85 61 41 25 13 5 c/3 144 112 84 60 40 24 12 4 b/3 17 15 13 11 9 7 5 3 a/3 此表与c-b=1时相同，为什么？ 本原勾股数 满足(a,b,c)=1的勾股数(a,b,c)称为本原勾股数。 由前面讨论知道，c-b=3时的勾股数不是本原勾股数，c-b=1时的勾股数是本原勾股数，c-b=2时如何呢？下面通过程序找出200以内的所有本原勾股数。 l=0; m=[]; for a=1:200 for c=a+1:100 b=sqrt(c^2-a^2); if (b==floor(b)) (gcd(gcd(a,b),c)==1) (ba) l=l+1; m(:,l)=[a,b,c,c-b]; end end end m 结果中，最后一行显示的是c-b，我们看看有无规律？ 这个规律便是：若{a,b,c}为本原勾股数,则有c-b为完全平方数或完全平方数的两倍. 6.Pell方程 ( D为正整数 ) 下面的两个pell方程，我们编程求1000以内的解 function m=pell(D,d,n)