实验三 Mathematica的综合程序设计

实验三:Mathematica的综合程序设计

1.验证角谷猜想.如果一个自然数,是偶数,就把它除以2,如果是奇数,就把它乘以3加1, 这样下去,经过若干次之后,必定得到数字1.试编写程序验证这一猜想.

In[ ]:= n=Input["n="];Print["n=",n];k=0;

While[n!=1,k++;If[EvenQ[n],n=n/2,n=3n+1]];

Print["k=",k]

Out[ ]= n= 23

k= 15

2.求一个整数的所有因子

方法一:因子竖着放的

In[ ]:= n=Input["n="];Print["n=",n];

For[i=1,i<=n,i++,If[Mod[n,i]==0,Print[i]]]

In[ ]:= n= 12

Out[ ]= 1

2

3

4

6

12

方法二:借助于列表,横着排放

In[ ]:= n=Input["n="];Print["n=",n];t={};

For[i=1,i<=n,i++,If[Mod[n,i]==0,AppendTo[t,i]]];t

In[ ]:= n= 12

Out[ ]= {1,2,3,4,6,12}

3.验证Nn=p1*p2*...*pn+1是不是素数

此题,n不宜过大,以免溢出.

In[ ]:= For[k=1;n=1,n<=30,n++,

k=k*Prime[n];

If[PrimeQ[k+1],Print["N[",n,"]=",k+1," is a Prime"],

Print["N[",n,"] =",k+1," is not a Prime"]]]

4.Fermat数是不是素数,欧拉证明了n=5不是.试试看.

n

In[ ]:= f[n_] := 221

In[ ]:= f[{1,2,3,4,5,6,7,8}]

Out[ ]= {5,17,257,65537,4294967297,18446744073709551617,340282366920938463463374 607431768211457,1157920892373161954235709850086879078532699846656405640394575840 07913129639937 }

In[ ]:= PrimeQ[%]

Out[ ]= {True,True,True,True,False,False,False,False}

5.求所有完全数

方法一

In[ ]:= For[n=2,n<=1000,n++,t={};

For[i=1,i<=n,i++,If[Mod[n,i]==0,AppendTo[t,i]]];t;

If[(Plus@@t)==2n,Print[n," yes"]]]

Out[ ]= 6 yes

28 yes

496 yes

方法二

In[ ]:= For[n=2,n<=1000,n++,s=0;

For[i=1,i<=n-1,i++,If[Mod[n,i]==0,s=s+i]]; If[s== n,Print[n," yes"]]] Out[ ]= 6 yes

28 yes

496 yes

In[ ]:= For[n=2,n<=1000,n++,s=0;

For[i=1,i<=n/2,i++,If[Mod[n,i]==0,s=s+i]]; If[s== n,Print[n," yes"]]] Out[ ]= 6 yes

28 yes

496 yes

我们来观察一下时间的长短

In[ ]:= For[n=2,n<=1000,n++,s=0;

For[i=1,i<=n-1,i++,If[Mod[n,i]==0,s=s+i]];

If[s==n,Print[n," yes"]]]//Timing

Out[ ]= 6 yes

28 yes

496 yes

{2.023 Second,Null}

In[ ]:= For[n=2,n<=1000,n++,s=0;

For[i=1,i<=n/2,i++,If[Mod[n,i]==0,s=s+i]];

If[s==n,Print[n," yes"]]]//Timing

Out[ ]= 6 yes

28 yes

496 yes

{1.382 Second,Null}

6.利用牛顿迭代法求下列方程的解

(1) x^3 - 2 x - 5 = 0, x0 = 2

In[ ]:= f[x_] := x^3 - 2 x - 5

x0 = 2.0;

For[i = 1, i<= 5, i++,x1 = x0 - f[x0] / f'[x0]; Print[i, " ", x1]; x0 = x1]

Out[ ]= 1 2.1

2 2.09457

3 2.09455

4 2.09455

5 2.09455

(2) x^41 + x ^3 - 1 = 0, x0 = 1

In[ ]:= f[x_] := x^41 + x ^3 - 1

x0 = 1.0;

For[i = 1, i<= 8, i++, x1 = x0 - f[x0] /f'[x0]; Print[i, " ", x1]; x0 = x1]

Out[ ]= 1 0.977273

2 0.960461

3 0.953391

4 0.952496

5 0.952484

6 0.952484

7 0.952484

8 0.952484

(3) x^3 - x - 1 = 0, x0 = 1.5

In[ ]:= f[x_] := x^3 - x - 1

x0 = 1.5;

For[i = 1, i<= 5,i++,x1 = x0 - f[x0] /f' [x0]; Print[i, " ", x1]; x0 = x1]

Out[ ]= 1 1.34783

2 1.3252

3 1.32472

4 1.32472

5 1.32472

(4) x^3 - x^2 - 1 = 0, x0 = 1.5

In[ ]:= f[x_] := x^3 -x^2 - 1

x0 = 1.5;

For[i = 1, i<=5,i++,x1 = x0 - f[x0] /f' [x0]; Print[i, " ", x1]; x0 = x1]

Out[ ]= 1 1.46667