数学实验课程实验报告

《数学实验》实验报告
学生姓名
学号
院系
专业
任课教师
二Ｏ一五年12 月9 日
南京信息工程大学实验（实习）报告实验课程实验名称第一次实验实验日期 2015-9-16 指导老师
专业年级姓名学号得分
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
实验目的：熟悉Mathematica软件包的使用。

实验内容：
1、用两种方式编写如下自定义函数，求在x=-2.0,x=1.0,x=5.0处的函数值，并画出函数
x在区间[-10，10]上的图像
代码如下：
f1=Plot[E^x*Sin[x],{x,-10,0}];
f2=Plot[Cos[x],{x,0,E}];
f3=Plot[Cos[x]*Sin[x],{x,-E,10}];
Show[f1,f2,f3];
以及：
f[x_/;x<0]:=E^x*Sin[x]
f[x_/;x>0&&x<E]:=Cos[x]
f[x_/;x>E]:=Cos[x]*Sin[x]
Plot[f[x],{x,-10,10}]
图像如下：
三条求值语句为：
f[-2.0]
f[1.0]
f[5.0]
函数值输出分别为：
-0.12306
0.540302
-0.272011
2、 分别用Plot3D, ParametricPlot3D 函数画出12
2
2
=++z y x （10,10<<<<y x ）
的图像。

1、 语句：
2、 图像：
3、 语句：
ParametricPlot3D[{Sin[u]*Cos[v],Sin[u]*Sin[v],Cos[u]},{u,0,Pi/2},{v,0,Pi/2}] 4、 图像：
Plot3D 1x ^2y^2,x,0,1,y,0,1
3、用Mathematica实现一个四人追逐问题，给出结果并划出追逐路线（如下图）。

语句：
v=1;t=18;dt=0.02;n=t/dt;
T={{{0,10}},{{10,10}},{{10,0}},{{0,0}}};
d=Sqrt[(x2-x1)^2+(y2-y1)^2];
For[j=1,j≤n,j++,For[i=1,i≤4,i++,x1=T[[i,j,1]];y1=T[[i,j,2]];
If[i≠4,x2=T[[i+1,j,1]];
y2=T[[i+1,j,2]],x2=T[[1,j,1]];
y2=T[[1,j,2]]];
x1=x1+v*dt*(x2-x1)/d;
y1=y1+v*dt*(y2-y1)/d;
T[[i]]=Append[T[[i]],{x1,y1}]]];
P=Graphics[{Line[T[[1]]],Line[T[[2]]],Line[T[[3]]],Line[T[[4]]],Line[{{0,10},{10,10},{10,0 },{0,0},{0,10}}]}];
Show[P,AspectRatio→1];
图像：
实验要求：
撰写实验报告
写出试验过程中所使用的Mathematica程序或语句和计算结果
南京信息工程大学实验（实习）报告实验课程数学实验实验名称第二次实验实验日期 2015-9-16 指导老师
专业年级姓名学号得分
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
实验目的：
练习π的求解方法。

实验内容：
4、用反正切函数的幂级数展开式结合有关公式求π，若要精确到以40位、50位数字，
试比较简单公式和Machin公式所用的项数。

（1）真实值：N[Pi,50]
（2）Arctan幂级数展开法:
40位：
k=100000;
S1=N[4*Sum[(-1)^(n-1)/(2n-1),{n,1,k}],40]
50位：
k=100000;
S1=N[4*Sum[(-1)^(n-1)/(2n-1),{n,1,k}],50]
（2）简单公式（有效位数为40）：
k = 10;
S = N[4*Sum[(-1)^(n - 1)*(1/2)^(2n - 1)/(2n - 1) + (-1)^(n - 1)*(1/3)^(2n - 1)/(2n - 1), {n, 1.k}], 40]
（3）简单公式（有效位数为50）：
k = 10;
S = N[4*Sum[(-1)^(n - 1)*(1/2)^(2n - 1)/(2n - 1) + (-1)^(n - 1)*(1/3)^(2n - 1)/(2n - 1), {n, 1.k}], 50]
（4）Machin公式（有效位数为40）：
k = 10;
S = N[4*Sum[4*(-1)^(n - 1)*(1/5)^(2n - 1)/(2n - 1) - (-1)^(n - 1)*(1/239)^(2n - 1)/(2n - 1), {n, 1.k}], 40]
（4）Machin公式（有效位数为50）：
k=10;
S=N[4*Sum[4*(-1)^(n-1)*(1/5)^(2n-1)/(2n-1)-(-1)^(n-1)*(1/239)^(2n-1)/(2n-1),{n,1.k}],50] 运行结果：
1.（1）3.1415926535897932384626433832795028841971693993751
（2）3.141582653589793488462643352029502893728
（3）3.1415826535897934884626433520295028937284193939649
（4）3.141592579606351210965104011041479112904
（5）3.1415925796063512109651040110414791129038298707062
（6）3.141592653589791696917279619620105448141
（7）3.1415926535897916969172796196201054481406519829326
综合结果可知，Machin公式所得结果比简单公式所得结果收敛的速度快。

5、用数值积分计算π，分别给出用梯形法和Simpson法精确到10位数字、用Simpson
法精确到15位数字时所用的项数n及的 近似值
梯形法：
n=5000;
y[x_]:=4/(1+x*x);
s1=(Sum[y[k/n],{k,1,n-1}]+y[0]+y[1]/2)/n;
N[s1,10]输出：3.141992647
Simpson法（精确到10位数字）：
n=5000;
y[x_]:=4/(1+x*x);
s2=(y[0]+y[1]+2*Sum[y[k/n],{k,1,n-1}]+4*Sum[y[(k-1/2)/n],{k,1,n}])/(6*n);
N[s2,10]
输出：3.141592654
Simpson法（精确到15位数字）：
n=5000;
y[x_]:=4/(1+x*x);
s2=(y[0]+y[1]+2*Sum[y[k/n],{k,1,n-1}]+4*Sum[y[(k-1/2)/n],{k,1,n}])/(6*n);
N[s2,15]
输出：3.14159265358979
6、用计算机模拟Buffon实验，给出n=1,000、10,000、1,000,000时的模拟结果。

(1) (n=1000时)
n=1000;a=20;l=10;
S4=Block[{i,m=0},
For[i=n,i>0,i--,m=m+If[Random[]*a/2<=l/2*Sin[Random[]*Pi/2],1,0]];N[(2*l*n)/(a*m),10]] 输出：3.105590062
(2)（n=10000时）
n=10000;a=20;l=10;
S4=Block[{i,m=0},
For[i=n,i>0,i--,m=m+If[Random[]*a/2<=l/2*Sin[Random[]*Pi/2],1,0]];N[(2*l*n)/( a*m),10]]
输出：
3.185727939
(3)（n=1,000,000时）
n=100000;a=20;l=10;
S4=Block[{i,m=0},
For[i=n,i>0,i--,m=m+If[Random[]*a/2<=l/2*Sin[Random[]*Pi/2],1,0]];N[(2*l*n)/( a*m),10]]
输出：
3.150697880
实验要求：
撰写实验报告
写出试验过程中所使用的Mathematica程序或语句和计算结果
南京信息工程大学实验（实习）报告实验课程数学实验实验名称第三次实验实验日期 2015-10-21 指导老师
专业年级姓名学号得分
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
实验目的：
熟悉差分方程的求解，以及相关金融问题的数学建模方法。

实验内容：
1、
试根据以上年率表，计算出每万元1~10年的月还款表。

程序如下：
1、五年以内：
For[a=10000;;r=0.0477/12;k=1,k≤5,k++,Print[k];
b=(a*r*((1+r)^(k*12)))/((1+r)^(k*12)-1);
Print[b//N]]
2、五年以上：
For[a=10000;;r=0.0504/12;k=6,k≤10,k++,Print[k];b=(a*r*((1+r)^(k*12)))/((1
+r)^(k*12)-1);Print[b//N]]
输出结果：
1、五年以内：
1 855.021
2 437.685
3 298.677
4 229.253
5 187.66
2、五年以上：
6 161.235
7 141.527
8 126.79
9 115.366
10 106.261
2、小李夫妇曾经准备申请商业贷款10万元用于购置住房，每月还款880.66元，25年
还清。

房产商介绍的一家金融机构提出：贷款10万元，每半月还款440.33元, 22年还清, 不过由于中介费手续费等原因，贷款时要预付4000元。

小李考虑，虽然预付费用不少，可是减少三年还款期意味着减少还款近3万2千元，
而每月多跑一趟，那不算什么.这机构的条件似乎还是蛮优惠的。

试通过计算两种贷款的利率水平，比较那种贷款更优惠。

程序如下：
1、商业贷款：
a=100000;k=25;b=880.66;r=.
FindRoot[(a*(r/12)*((1+(r/12))^(k*12)))/((1+(r/12))^(k*12)-1)-b,{r,0.01}]
2、金融机构贷款：
a=100000-4000;k=2*22;b=440.33;r=.
FindRoot[(a*(r/24)*((1+(r/24))^(k*12)))/((1+(r/24))^(k*12)-1)-b,{r,0.01}]
输出结果：
1、商业贷款：
{Null (r→0.0960003)}
2、金融机构贷款：
{Null (r→0.0969953)}
综合得出，由于金融机构贷款利率高于商业贷款利率，所以商业利率贷款更优惠。

3、试比较两种提前还款方式的优劣（附加）
所谓提前还贷是指借款人在保证按月按额偿还个人住房贷款本息的基础上，提前偿
还部分或全部购房借款的一种经济行为。

每次提前还款后，相应冲减余贷款本金。

银行根据尚未归还的贷款本金重新计算借款人的月均还款额，直至贷款本息全部还
清。

重新计算月还款金额有两种方式：
A、提前还款额冲抵最后月份的本金，每月的还款额度不变，还款时间缩短；
B、提前还款额冲抵本金后，将剩余的贷款重新计算月还款额减少，还款时间不变。

例如，谢先生申请公积金贷款30万元，贷款期限为20年，在正常按月还了5年贷
款后，谢先生决定提前还5万元本金，然后再继续按月还款。

试比较两种提前还款方式的优劣？‘
程序如下：
A方案：
a=300000;r=0.0504/12;k=20;
b=(a*r*((1+r)^(k*12)))/((1+r)^(k*12)-1);n=b*15*12-50000;
Print[n//N] ;
a=300000;r=0.0504/12;k=20;
b=(a*r*((1+r)^(k*12)))/((1+r)^(k*12)-1);n=50000/b;
Print[n//N]
B方案：
a=300000;r=0.0504/12;k=15;
b=15*12*((a*15/20-50000)*r*((1+r)^(k*12)))/((1+r)^(k*12)-1)+50000;
Print[b//N];
输出结果：
A方案：
307570.
25.1699
B方案：
299757.
综合得出，A方案中，还款额大于B方案，但是能提前两年半把余款还
清。

B方案还是继续还款15年，但是还款总额较少。

实验要求：
撰写实验报告
写出试验过程中所使用的Mathematica程序或语句和计算结果
南京信息工程大学 实验（实习）报告
实验课程 数学实验 实验名称 第四次实验报告 实验日期 2015-10-28 指导老师 专业 年级 姓名 学号 得分
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
实验目的：
熟悉迭代法的基本概念，并用迭代法求解方程、方程组的根。

实验内容：
1、 自己构造2种不同的迭代格式求32的近似值，并比较收敛速度。

2、 第一题： 方法一：（牛顿迭代法） 语句：
f[x_]:=x^3-2;g[x_]:=3x^2 ;
x0=-2.;esp=10^(-6);For[i=1,i ≤11,i++,x1=x0-f[x0]/g[x0];If[Abs[x1-x0]>esp,x0=x1,Break[]];Print[x1];]输出结果：
-1.16667 -0.287982 7.84658 5.24188 3.51885 2.39974 1.171559 1.37023 1.26856 1.25998 1.25992 方法二：（弦位法） 语句：f[x_]:=x^3-2 ; Plot[f[x],{x,-2,2}] ;
FindRoot[f[x],{x,0,1}] [f[x],{x,0,1}] 输出结果：
{x →1.25992} 收敛速度相同。

3、 画出2cos sin x x x =的图像，并利用牛顿迭代法求出该方程的所有根。

语句：
f[x_]:=Sin[x]*Cos[x]-x^2;D[f[x],x] 输出：
2x Cos x
2
Sin x
2
f[x_]:=Sin[x]*Cos[x]-x^2;
g[x_]:=-2 x+Cos[x]^2-Sin[x]^2;
Plot[f[x],{x,-2,2}] ;
x0=0.4;esp=10^(-10);
For[i=1,i ≤10,i++,x1=x0-f[x0]/g[x0];If[Abs[x1-x0]>esp,x0=x1,Break[]];Print[x1];]
输出结果：
2.32344 1.07198 0.800582 0.71406
0.71406 0.702425 0.702207 0.702207
4、 对方程组b Ax =，设A 的对角元素0≠ii a ，
令),,,(2211nn a a a Diag D =为对角阵，
将方程组改写成b x A D Dx +-=)(，或b D x A D I x 1
1)(--+-=
用这种迭代格式求解方程组b Ax =，其中 ⎪⎪⎪⎭
⎫ ⎝⎛--=21111
1112A ，b=0 并将结果与f Mx x +=迭代格式的结果进行比较。

语句：
Jacob 迭代格式：
SeideIterate[a_,b_List,x0_List,n_Integer]:=Module[{ad=Length[a],i,j,k,var=x0},For[i=1,i<=ad,i++,If[a[[i,i]]==0,Print["a[",i,",",i,"]=0."];
Abort[]]];
For[i=1,i<=n,i++,Print[var];
For[j=1,j<=ad,j++,var[[j]]=N[(b[[j]]-Sum[a[[j,k]]*var[[k]],{k,ad}])/a[[j,j]]+var[[j]],20] ];
]
a={{2,-1,1},{1,1,1},{1,1,-2}};b={0,0,0};x0={1,1,1};
SeideIterate[a,b,x0,20];
输出结果：
{1,1,1} {0.,-1.,-0.5} {-0.25,0.75,0.25} {0.25,-0.5,-0.125}
{-0.1875,0.3125,0.0625} {0.125,-0.1875,-0.03125} {-0.078125,0.109375,0.015625} {0.046875,-0.0625,-0.0078125}{-0.0273438,0.0351563,0.00390625}
{0.015625,-0.0195313,-0.00195313}
{-0.00878906,0.0107422,0.000976563}
{0.00488281,-0.00585938,-0.000488281}
{-0.00268555,0.00317383,0.000244141}
{0.00146484,-0.00170898,-0.00012207}
{-0.000793457,0.000915527,0.0000610352}
{0.000427246,-0.000488281,-0.0000305176}
{-0.000228882,0.000259399,0.0000152588}
{0.00012207,-0.000137329,-7.62939×10-6}
{-0.0000648499,0.0000724792,3.8147×10-6}
{0.0000343323,-0.000038147,-1.90735×10-6}
Seidel迭代格式：
语句：
LSIterate[m_,f_List,f0_List,n_Integer]:=
Module[{i,var=f0,t=Table[{},{i,n}]},
For[i=1,i£n,i++,t[[i]]=var;var=m.var+f];
t]
m={{0.33,0.11,0.22},{-0.33,0.56,0.11},{0,0.33,-0.33}};
f={1,1,1};f0={0,0,0};
LSIterate[m,f,f0,25]
输出结果：
{{0,0,0},{1.,1.,1.},{1.66,1.34,1.},{1.9152,1.3126,1.1122},
{2.02109,1.22538,1.06613},{2.0363,1.13653,1.05255},
{2.02856,1.08026,1.02771},{2.01435,1.04857,1.01734},
{2.00389,1.03437,1.01031},{1.99733,1.0291,1.00794},
{1.99407,1.02805,1.00698},{1.99266,1.02843,1.00695},
{1.99224,1.02911,1.00709},{1.9922,1.02964,1.00727},
{1.99229,1.02997,1.00738},{1.99238,1.03014,1.00745},
{1.99244,1.03022,1.00749},{1.99248,1.03024,1.0075},
{1.99249,1.03024,1.0075},{1.9925,1.03024,1.0075},
{1.9925,1.03023,1.0075},{1.9925,1.03023,1.0075},
{1.9925,1.03023,1.0075},{1.9925,1.03023,1.0075}, {1.9925,1.03023,1.0075}}
实验要求：
撰写实验报告
写出试验过程中所使用的Mathematica程序或语句和计算结果
南京信息工程大学 实验（实习）报告 实验课程 数学实验 实验名称 第五次实验 实验日期 2015-11-11 指导老师 专业 年级 姓名 学号 得分 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
实验目的：
了解有关分形和混沌的基本理论，能够用Mathematica 软件绘制出一些简单的分形和混沌图形。

实验内容：
1、 用Mathematica 软件绘制一个分形的图形，图形类别自选。

Koch 雪花曲线
程序如下：
redokoch[ptlist_List] := Block[{tmp = {}, i, pnum = Length[ptlist]},
For[i =1, i < pnum, i++,
tmp = Join[tmp, {ptlist[[i]], ptlist[[i]]*2/3 + ptlist[[i + 1]]/3,
(ptlist[[i]] + ptlist[[i + 1]])/2 + {
ptlist[[i]][[2]] - ptlist[[i + 1]][[2]],
ptlist[[i + 1]][[1]] - ptlist[[i]][[1]]}*Sqrt[3]/6,
ptlist[[i]]/3 + ptlist[[i + 1]]*2/3, ptlist[[i + 1]]}]];
tmp] ; lnko01 = {{0, 0}, {1, 0}}; Show[
Graphics[Line[Nest[redokoch, lnko01, 5]], AspectRatio -> Sqrt[3]/6]]
图像如下：
2、 令⎩⎨⎧-==+==1
)(5.01)(5.02211sz z p sz z p ωω，其中I s 5.05.0+=，绘制出相应的IFS 吸引子图形，并
取不同的s ，观察图形的变化。

程序如下：
s=0.5+0.5*I;
p1=0.5;f1[z_]:=s*z+1;
p2=0.5 ;f2[z_]:=s*z-1;
f[z_]:=Block[{tmp},tmp=Random[];
Which[tmp<p1,f1[z],tmp<1,f2[z]]];
Array[mu,{150,150}];
showIFS[z0_,shrage_List,divi_List,nmax_]:=
Block[{i,j,z=z0,a=divi[[1]],b=divi[[2]],temp1,temp2,mumax=0},
For[i=a,i>=1,i--,For[j=b,j>=1,j--,mu[i,j]=0]];
For[i=nmax,i>=1,i--,temp1=
Floor[a*(Re[z]-shrage[[1]][[1]])/(shrage[[2]][[1]]-shrage[[1]][[1]])]+1;temp2=Floor [b*(Im[z]-shrage[[1]][[2]])/(shrage[[2]][[2]]-shrage[[1]][[2]])]+1;mu[temp1,temp2] ++;
z=f[z];];For[i=a,i>=1,i--,For[j=b,j>=1,j--,mumax=Max[mumax,mu[i,j]]]];
mu1=Table[GrayLevel[1-N[mu[j,i]]/mumax],{i,a},{j,b}];
Show[Graphics[RasterArray[mu1]]]] ;
showIFS[0+I 0,{{-0.1,-0.1},{1.1,1.1}},{150,150},10000]
图像如下：
3、用Mathematica软件绘制一个混沌的图形，图形类别自选.
用二次迭代序列迭代函数f(x)=ax(1-x)
程序如下：
IterGeo[u_, x0_] := Module[{p1, p2, i,
pointlist = {},
var = x0, fvar = u*x0*(1 - x0)}, p1 = Plot[{u*x*(1 - x), x}, {x, 0, 1},
DisplayFunction -> Identity]; For[i = 1, i < 20,
i++, AppendTo[pointlist, {var, fvar}]; AppendTo[
pointlist, {fvar, fvar}]; var = fvar; fvar = u*var*(1 - var)]; p2 =
ListPlot[pointlist,
PlotJoined -> True, DisplayFunction -> Identity];
Show[{p1, p2}, DisplayFunction -> $DisplayFunction]]; IterGeo[3.6, 0.8]
图像：
4、谈谈你所认识的分形和混沌。

答：分形：具有无限嵌套层次的精细结构，且在不同尺度下保持相似属性。

混沌：对初值敏感，而且不是随机的。

实验要求：
撰写实验报告
写出试验过程中所使用的Mathematica程序或语句和计算结果。