数学实验报告基于matlab.docx

数学实验报告基于matlab.docx

《数学实验报告基于matlab.docx》由会员分享，可在线阅读，更多相关《数学实验报告基于matlab.docx（38页珍藏版）》请在冰豆网上搜索。

数学实验报告基于matlab

西安交通大学

实验报告

课程：

高数实验

班级：

自动化32

姓名：

张文杰

学号：

2130504049

日期：

2014年4月20日

1、

计算值其它公式：

请按照上述公式分别编程计算并与

（1）

就精度与叠加次数进行比较，能得出什么结论？

？

解：

（1）按照以下的程序得到n取不同的值

（1）式计算的

的近似值

formatlong

x=1:

6;%首先计算原函数的由10的1到6次方展开时的值

forn=10.^x%利用循环一次性将公式展开为不同的次数

s=0;

digits（25）;

fork=1:

n

s=s+4*（-1）^（k+1）/（2*k-1）;%计算出的近似值

end

vpa（s,25）%定义精度为25位

end

（2）分别计算a,b,c,d四个公式的

的近似值

x=1:

6;

form=10.^x%计算上面四个公式的不同展开次数的值

fa=0;fb=0;fc=0;fd=0%分别表示四个等式的值

digits（25）;

fork=1:

m

fa=fa+1/（2*k-1）^2;

fb=fb+（-1）^（k-1）/k^2;

fc=fc+（-1）^（k-1）/（2*k-1）^3;

fd=fd+sin（2*k-1）/（2*k-1）^3;

end

fa=vpa（sqrt（8*fa）,25）

fb=vpa（sqrt（12*fb）,25）

fc=vpa（（32*fc）^（1/3）,25）

fd=vpa（（1+sqrt（1+32*fd））/2,25）%

end%

（3）将体现四个公式收敛情况的图形画到一个直角坐标系上

pfa=[];pfb=[];pfc=[];pfd=[];ps=[];%绘制出体现收敛速度的图像

t=100;%选取的项数为150

step=5;%步长选为5

forj=10:

step:

t%

fa=0;fb=0;fc=0;fd=0;s=0;

fork=1:

j

fa=fa+1/（2*k-1）^2;

fb=fb+（-1）^（k-1）/k^2;

fc=fc+（-1）^（k-1）/（2*k-1）^3;

fd=fd+sin（2*k-1）/（2*k-1）^3;

s=s+4*（-1）^（k+1）/（2*k-1）;

end

pfa=[pfa,sqrt（8*fa）];

pfb=[pfb,sqrt（12*fb）];

pfc=[pfc,（32*fc）^（1/3）];

pfd=[pfd,（1+sqrt（1+32*fd））/2];%

ps=[ps,s];

end

j=10:

step:

t;

plot（j,pfa,'b',j,pfb,'g',j,pfc,'r',j,pfd,'y',j,ps,'k'）%绘制收敛的图像

程序的运行结果：

叠加次数

函数

结果

函数

结果

原函数的

不同展开

次数的值

ans=

.0418********

ans=

3.131********8553686593768

ans=

3.1405926538397941349956

ans=

3.141492653590034489496929

ans=

3.141582653589719775766298

ans=

3.141591653589774324473183

N=10

fa=

3.10962545798864775647985

fb=

3.132********9482278866735

fc=

3.141526087929505717255552

fd=

3.141538602696021698079676

N=100

fa=

3.138********7091242086008

fb=

3.141498114035649713571274

fc=

3.141592586052465385648702

fd=

3.141592599756393155985279

N=1000

fa=

3.141274327602740168430273

fb=

3.141591699614915800253812

fc=

3.141592653522246259001349

fd=

3.141592653631021470062024

N=10000

fa=

3.141560822439948719164704

fb=

3.141592644041437232260705

fc=

3.141592653589793238462643

fd=

3.141592653589710071315722

N=

100000

fa=

3.141589470489344115833319

fb=

3.141592653494265530156326

fc=

3.141592653589793238462643

fd=

3.141592653589793238462643

N=

1000000

fa=

3.141592335279969727679372

fb=

3.141592653588781924867135

fc=

3.141592653589793238462643

fd=

3.141592653589793238462643

结果分析：

当叠加次数N=100000时五个公式分别精确到第3、4、7、14、13位，由此可见当叠加次数相同的情况下，精确度的顺序由大到小为c>d>b>a>

（1）式，这也正好验证了幂级数展开式中级数越高，收敛的速度也就越快，结果也就越精确。

绘图运行结果：

结果分析：

已设定a\b\c\d分别为黑色、绿色、红色、黄色，由结果可以直观的看出c的收敛速度是最快的，a是最慢的与以上的数字实验结果一致。

实验体会：

本次实验没有按照书上的程序分别代入不同的n值来进行多次运行，而是采取了for循环，一次性算出六个不同展开次数的精确值，节约了大量的时间。

在将运行结果导入报告时，直接用excel将数据导入了表格（相当的方便啊，，以前都没发现，一个一个的输入）

过程总结：

做的时候有还好几处都是参数弄混了，也报不出来错，但就是运行不出来，那个汗哪！

！

！

！

！

！

！

！

！

！

！

！

！

为了避免麻烦，将所有的程序放在一起运行，一起出结果，所以结果很直观。

2、基于关系式

，利用蒙特卡洛方法近似计算

解：

运用蒙特卡洛方法编制程序如下：

n=input（'请输入产生点的个数：

'）

m=0;

fori=1:

n

a=rand（1,2）;

ifa

（1）^2+a

（2）^2<=1

m=m+1;

end

end

mypi=4*m/n

运行结果：

>>pi

请输入产生点的个数：

1000

n=1000

mypi=3.104000000000000

请输入点的个数

N=

Mypi=

1000

1000

3.104000000000000

5000

5000

3.132********0000

10000

10000

3.142400000000000

50000

50000

3.136********0000

100000

100000

3.146760000000000

结果分析：

通过依次增大输入点的个数观察结果发现，这种方法的精确度无法保证，点的个数越多，不一定就精确度越高，而且计算的速度超慢（试过1000000000）老半天没反应，只能弃疗。

3、设有一制作均匀的冰淇淋可以看做由圆锥面和球面围成，采用取随机数的方法，用蒙特卡洛法计算这个冰淇淋的体积。

解：

程序：

s=0;

n=input（'请输入产生点的个数：

'）;

fori=1:

n

a=rand（1,2）-1;

b=rand（1,1）*2;

ifsqrt（a

（1）^2+a

（2）^2）<=b

（1）&&b

（1）<=1+sqrt（1-a

（1）^2-a

（2）^2）

s=s+1;

end

end

v=8*s/n

运行结果：

请输入产生点的个数：

10000v=3.206400000000000

请输入产生点的个数：

10000v=3.139********0000

请输入产生点的个数：

10000v=3.133********0000

请输入产生点的个数：

10000v=3.160800000000000

结果分析：

算出来后不知对不对，就用理论算法做出来是

，之后就肆无忌惮的带值算结果真的好随机啊，以上是我带相同的点进去，差别是如此的大

总结：

对了，这次试验中间出了点小事故由于不清楚rand（1,2）只产生1——2之间的随机数，吃了不少苦头，检查了半天算出来老是零点几。

。

。

。

。

。

。

。

。

练习八

4、某医院每日至少需要护士人数如表所示

每班护士在值班开始时向病房报道，连续工作8小时，医院至少需要多少护士才能满足值班要求？

解：

（1）问题分析：

这是一个求最优解的问题，满足的两个大条件是每人工作8小时，且每个时间段都有人数限制，现在可以将所有24小时的班展现在同一时刻，每个时刻的人数均大于等于规定人数即可。

建立模型：

设第i个班次开始上班的护士数为

，建立以下模型

编写程序：

c=[1,1,1,1,1,1];

a=[-1,0,0,0,0,-1;-1,-1,0,0,0,0;

0,-1,-1,0,0,0;0,0,-1,-1,0,0;

0,0,0,-1,-1,0;0,0,0,0,-1,-1];

b=[-60;-70;-60;-50;-20;-30];

x=linprog（c,a,b）

miny=c*x

程序运行的结果：

x=

42.0000

28.0000

35.0000

15.0000

10.0000

20.0000

miny=150.0000

结果分析：

由此可以得出护士的总人数应至少为150人。

5、某工厂利用甲，乙两种原料生产A1，A2，A3三种产品。

每月可供应的原料数量（单位：

t），每万件产品所需各种原料的数量及每万件产品的价格如表所示：

试制定每月的最优生产计划，使得总收益最大。

解：

与上题差不多，就是简单地线性规划求解的问题

建立模型：

设三种产品的数量分别为

，

Maxy=

编写下列程序：

c=[-12;-5;-4];

a=[4,3,1;2,6,3];

b=[180;200];

aeq=[];

beq=[];

vlb=[0;0;0];

vub=[];

[x,fval]=linprog（c,a,b,aeq,beq,vlb,vub）

maxy=-fval

运行结果：

x=

34.0000

0.0000

44.0000

fval=

-584.0000

maxy=

584.0000

结果分析：

当生产34万件A1产品，，44万件A2产品是总收益最大为584万元。

总结：

这个题本来挺简单地，也明白linrog函数只能求最小值，做的时候只知道把最后结果乘以负号但是C矩阵给忘了，导致最后结果出现了复数，不过很快就纠正过来了。

6、一公司拟在某市东，南，西三区建立门市部，拟议中有7个位置

可供选择，规定东区在

三点中至多选两个。

西区在

两点中至少选一个，南区在

两个点中至少选一个，并知道如选用

点，则投资为

元，估计每年可获利润为

元，但投资总额不超过B元，问应选择哪几个点可使年利润最大？

解：

假设当选中

时对应的投资额为30,20,40,10,50,60,40.对应的可获利润为，110,140,130,150,190,130,170，总投资不超过150，所以可以建立以下的数学模

编制以下程序：

c=[-110,-140,-130,-150,-190,-130,-170];

a=[30,20,40,10,50,60,40;

1,1,1,0,0,0,0;

0,0,0,-1,-1,0,0;

0,0,0,0,0,-1,-1];

b=[150;2;-1;-1];

aeq=[];

beq=[];

vlb=[0;0;0;0;0;0;0];

vub=[1;1;1;1;1;1;1];

x=linprog（c,a,b,aeq,beq,vlb,vub）

shouyi=d*x

运行结果

x=1.0000

1.0000

0.0000

1.0000

1.0000

0.0000

1.0000

shouyi=

760.0000

结果分析：

本题由自己设置数据，做出来的结果貌似还不错，即选择

建立门市部所得的年利润最大。

总结：

通过自己设定数据，让我对于在实际生活中怎样把数据转化为可供计算的数学问题，进一步提高了解决实际问题的能力，还有就是这次程序一次成功，说明是可以熟能生巧的。

7、某单位有三百万元可用于投资，共有六个项目可供选择，其投资额分别为40，60，80，50，90，70（万元），预计三年可获利润分别为10，12，15，11，16，13（万元），试确定一种投资可使三年获得的利润最大？

解：

问题分析：

投资方式的挑选要满足在总投资额小于300万元的前提下，取得最大收益值，即可视为一个0-1规划问题。

1，第j种方式被选择

令

0，第j种方式未被选

该问题的数学模型是：

maxZ=

实验程序：

k=1;

zl=[10,12,15,11,16,13];

sy=[40,60,80,50,90,70];

forg1=0:

1

forg2=0:

1

forg3=0:

1

forg4=0:

1

forg5=0:

1

forg6=0:

1

fa=[g1,g2,g3,g4,g5,g6];

ifsy*fa'<=300

fazl（k,1:

6）=fa;

k=k+1;

end

end

end

end

end

end

end

falr=fajz*zl';

[f,i]=max（falr）

xuanfa=fazl（i,1:

6）

fatz=sy*xuanfa'

运行结果：

f=

61

i=

58

xuanfa=

111101

结果分析：

利用0—1规划得，此公司可选择除第五种外的投资类型从而获得最大利润，

最大利润为61万元。

8、某企业在两个相互分离的市场上出售同一产品，两个市场的需求函数分别为其中分别表示该产品两个市场的价格（单位：

万元/吨），分别表示该产品在两个市场的销售量（单位：

吨）.该企业生产这种产品总成本函数为其中q表示该产品在两个市场的销售总量，即,在产销平衡的状态下：

（1）如果该企业实行价格差别策略，试确定两个市场上该产品的销售量和最优价格，使该企业获得最大利润；

（2）如果该企业实行价格无差别策略，试确定两个市场上该产品的销售量和最优价格，使该企业获得最大利润，并比较两种价格策略下总利润的大小。

9、表给出3只股票49个周末的收盘价，

10、Feigenbaum在做实验的时候，对超越函数y=r*sin（pi*x）（r为非负实数）进行了分叉与混沌的研究，试用迭代格式x（i）=r*sin（pi*x（i-1））,做出相应的Feigenbaum图。

思路：

用axis，grid绘制网状坐标轴，以0.005对r进行微小改变，用fprintf指令绘图，并holdon观察图形即可。

程序：

holdon

axis（[0.5,1,0,1]）;grid

forr=0.5:

0.005:

1

x=[0.1];

fori=2:

150

x（i）=r*sin（pi*x（i-1））;

end

pause（0.001）

fprintf（'r=%.3f\n',r）

fori=101:

150

plot（r,x（i）,'k.'）;

end

end

运行结果：

结果分析：

如图即为该函数所做出的图像

11、求非线性方程组[2x

（1）^2-x

（1）*x

（2）-5*x

（1）+1;x

（1）+3lgx

（1）-x

（2）^2]的解，初值如下：

（1）x0=[1.4,-1.5];

（2）x0=[3.7,2.7]

思路：

先在文件编辑器中写下函数文件，再再命令框中使用求解命令。

程序：

函数文件

functionf=group1（x）

f=[2*x

（1）^2-x

（1）*x

（2）-5*x

（1）+1;

x

（1）+3*log（x

（1））-x

（2）^2]

命令框中：

[x,fval]=fsolve（'group1',[1.4,-1.5]）

结果：

x=

1.3735-1.5250

fval=

1.0e-006*

0.0313

-0.2032

命令框中：

[x,fval]=fsolve（'group1',[3.7,2.7]）

结果：

x=

3.75682.7798

fval=

1.0e-012*

0.1066

-0.4512

12、有一形状较为复杂，但表面很光滑的曲面工件，通过科学手段，将其放置于某一空间坐标系下，测得曲面上若干个点的坐标如下：

Xy

-5

-4

-3

-2

-1

0

1

2

3

4

5

-5

13.6

-8.2

-14.8

-6.6

1.4

0

-3.8

1.4

13.6

16.8

0

-4

-8.2

-15.8

-7.9

2.2

3.8

0

0.6

7.3

10.1

0

-16.8

-3

-14.8

-7.9

2.5

5.8

2.3

0

2.7

5.1

0

-10.1

-13.7

-2

-6.6

2.2

5.9

3.0

-0.3

0

1.9

0

-5.1

-7.3

-1.4

-1

1.4

3.8

2.3

-0.3

-0.9

0

0

-1.7

-2.7

-0.6

3.8

0

0

0

0

0

0

0

0

0

0

0

0

1

-3.8

0.6

2.7

1.7

0

0

0.9

0.3

-2.3

-3.8

-1.4

2

1.4

7.3

5.1

0

-1.7

0

0.3

-3.1

-5.8

-2.2

6.6

3

13.6

10.1

0

-5.1

-2.7

0

-2.3

-5.8

-2.5

7.9

14.8

4

16.8

0

-10.1

-7.3

-0.6

0

-3.8

-2.2

7.9

15.8

8.2

5

0

16.3

-13.6

-1.4

3.8

0

-1.4

6.6

14.8

8.2

-13.6

程序如下：

x=-5:

1:

5;

y=-5:

1:

5;

[xb,yb]=meshgrid（x,y）;

zb（1,:

）=[13.6,-8.2,-14.8,-6.6,1.4,0,-3.8,1.4,13.6,16.8,0];

zb（2,:

）=[-8.2,-15.8,-7.9,2.2,3.8,0,0.6,7.3,10.1,0,-16.8];

zb（3,:

）=[-14.8,-7.,2.5,5.8,2.3,0,2.7,5.1,0,-10.1,-13.7];

zb（4,:

）=[-6.6,2.2,5.9,3.0,-0.3,0,1.9,0,-5.1,-7.3,-1.4];

zb（5,:

）=[1.4,3.8,2.3,-0.3,-0.9,0,0,-1.7,-2.7,-0.6,3.8];

zb（6,:

）=[0,0,0,0,0,0,0,0,0,0,0];

zb（7,:

）=[-3.8,0.6,2.7,1.7,0,0,0.9,0.3,-2.3,-3.8,-1.4];

zb（8,:

）=[1.4,7.3,5.1,0,-1.7,0,0.3,-3.1,-5.8,-2.2,6.6];

zb（9,:

）=[13.6,10.1,0,-5.1,-2.7,0,-2.3,-5.8,-2.5,7.9,14.8];

zb（10,:

）=[16.8,0,-10.1,-7.3,-0.6,0,-3.8,-2.2,7.9,15.8,8.2];

zb（11,:

）=[0,16.3,-13.6,-1.4,3.8,0,-1.4,6.6,14.8,8.2,-13.6];

mesh（xb,yb,zb）

xc=-5:

0.25:

5;

yc=-5:

0.25:

5;

[xcb,ycb]=meshgrid（xc,yc）;

zcb=interp2（xb,yb,zb,xcb,ycb）;

mesh（xcb,ycb,zcb）

（1）

（2）所有点的竖坐标：

zcb=

Columns1through6

13.6000

8.1500

2.7000

-2.7500

-8.2000

-9.8500

8.1500

3.5875

-0.9750

-5.5375

-10.1000

-10.8438

2.7000

-0.9750

-4.6500

-8.3250

-12.0000

-11.8375

-2.7500

-5.5375

-8.3250

-11.1125

-13.9000

-12.8313

-8.2000

-10.1000

-12.0000

-13.9000

-15.8000

-13.8250

-9.8500

-10.7875

-11.7250

-12.6625

-13.6000

-11.5250

-11.5000

-11.4750

-11.4500

-11.4250

-11.4000

-9.2250

-13.1500

-12.1625

-11.1750

-10.1875

-9.2000

-6.9250

-14.8000

-12.8500

-10.9000

-8.9500

-7.0000

-4.6250

-12.7500

-10.7375

-8.7250

-6.7125

-4.7000

-2.6875

-10.7000

-8.6250

-6.5500

-4.4750

-2.4000

-0.7500

-8.6500

-6.5125

-4.3750

-2.2375

-0.1000

1.1875

-6.6000

-4.4000

-2.2000

0.0000

2.2000

3.1250

-4.6000

-2.8000

-1.0000

0.8000

2.6000

3.2000

-2.6000

-1.2000

0.2000

1.6000

3.0000

3.2750

-0.6000

0.4000

1.4000

2.4000

3.4000

3.3500

1.4000

2.0000

2.6000

3.2000

3.8000

3.4250

1.0500

1.5000

1.9500

2.4000

2.8500

2.5687

0.7000

1.0000

1.3000

1.6000

1.9000

1.7125

0.3500

0.5000

0.6