中南大学数学建模实验报告

MATLAB数学建模实验报告

学院：材料科学与工程

专业班级：材料国际

姓名：

学号：

完成时间：2016年12月7日

目录

一、数学实验学习体会 (3)

二、实验一：MATLAB作图 (4)

实验目的 (4)

实验内容 (4)

三、实验二：线性规划 (7)

实验目的 (7)

实验内容 (7)

四、实验三：插值 (11)

实验目的 (11)

实验内容 (11)

五、实验四：拟合 (12)

实验目的 (12)

实验内容 (12)

六、实验五：MATLAB在材料力学里的应用 (14)

实验目的 (14)

实验内容 (15)

七、实验六：MATLAB创建2048小游戏 (19)

游戏规则 (20)

游戏代码及运行结果 (20)

八、心得与收获 (26)

一、数学实验学习体验

通过对《数学实验与建模》这门课程的学习，我初步掌握了一些建模思想、模型分析以及对于数学矩阵实验室（即：MATLAB软件）的使用。课程分为两个阶段，即八周的数学建模讲授、八周的数学实验。在这里，主要谈一谈运用MATLAB软件进行的数学实验给我带来的感受与收获。

通过学习，我们知道MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。正如这些强大的功能，我们才有必要认真学习并掌握这项技能，我的专业是材料科学与工程，矩阵运算和处理实验数据对于我的专业领域大有裨益，这也坚定了我熟练掌握MATLAB的决心。

我做的第一个实验是图形的绘制。这在Microsoft软件中也可以实现，而MATLAB给我带来的直观感受就是更加“强大、丰富、专业”，不仅包含了二维三维，甚至多维度空间图形也能表现出来。还可对坐标控制、图形修饰、窗口分割等操作，如果特殊需要时还可用polar得到极坐标图形，调用semilogx得到对数坐标函数等。三维图形有三维曲线、三维曲面，这种功能对求两个复杂三维立体图形的交线交面等很有帮助。在二维图形绘制时可以绘出条形图、杆图、饼图，当然也可以调用函数bar3、stem3、pie3、fill3等绘制三维图形。对三维图形可以进行精细处理，比如视点处理，色彩处理，还可以进行图形的裁剪，在实际生活中也很有用。

另外一个让我影响深刻的功能就是数据处理，对于材料科学的科研工作，往往需要在大量实验数据里找到一定规律，从而揭示一种材料性能的影响因素，实现对材料性能的调控。而从MATLAB中最初学习到的就是插值与拟合，种类丰富，处理也十分精确，还可以自定义插值、拟合函数，最后通过plot以图形的形式展现出来。对于数据规律性的探讨十分有帮助。

通过这么短时期的学习，是很难理解到MATLAB的精髓的，要想从使用到理解到熟练掌握还需要一个很长的学习探索过程，我相信，MATLAB软件不仅将对我的科研领域起到重要的作用，还将为我处理生活问题带来便捷。

二、实验一：MATLAB作图

1.实验目的：

了解MATLAB作图的基本内容

掌握MATLAB作图的集中基本方式，实现对数据进行可视化作图分析2.实验内容：

（1）例：在[0,2pi]用红线画sin x,用绿圈画cos x.

x=linspace(0,2*pi,30);

y=sin(x);

z=cos(x);

plot(x,y,'r',x,z, 'g0')

先建M文件myfun1.m：

function Y=myfun1(x)

Y=exp(2*x)+sin(3*x.^2)

再输入命令：

fplot('myfun1',[-1,2])

(3)例:画多条曲线观察函数Z=(X+Y)2 x=-3:0.1:3;y=1:0.1:5;

[X,Y]=meshgrid(x,y);

Z=(X+Y).^2;

plot3(X,Y,Z)

(4)例:画函数Z=(X+Y)2的图形.

x=-3:0.1:3;

y=1:0.1:5;

[X,Y]=meshgrid(x,y);

Z=(X+Y).^2;

surf(X,Y,Z)

shading flat

(5)例:在区间[0,2π]画sin(x)的图形，并加注图例“自变量X”、“函数Y”、“示

意图”, 并加格栅.

x=linspace(0,2*pi,30);

y=sin(x);

plot(x,y)

xlabel('自变量X')

ylabel('函数Y')

title('示意图')

grid on

(6)例:山峰的三维和二维等值线图

[x,y,z]=peaks;

subplot(1,2,1)

contour3(x,y,z,16,'s')

grid, xlabel('x-axis'),ylabel('y-axis')

zlabel('z-axis')

title('contour3 of peaks');

subplot(1,2,2)

contour(x,y,z,16,'s')

grid, xlabel('x-axis'), ylabel('y-axis')

title('contour of peaks');

三、实验二：线性规划

1.实验目的：

了解线性规划的基本内容.

掌握用数学软件包求解线性规划问题

2.实验内容：

(1)例： max 6543216.064.072.032.028.04.0x x x x x x z +++++= 123456s.t.0.010.010.010.030.030.03850x x x x x x +++++≤ 70005.002.041≤+x x 10005.002.052≤+x x 90008.003.063≤+x x

1,2,,6j x j ≥=

解：编写M 文件xxgh1.m 如下：

c=[-0.4 -0.28 -0.32 -0.72 -0.64 -0.6];

A=[0.01 0.01 0.01 0.03 0.03 0.03;0.02 0 0 0.05 0 0;0 0.02 0 0 0.05 0;0 0 0.03 0 0

0.08];

b=[850;700;100;900]; Aeq=[]; beq=[];

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

[x,fval]=linprog(c,A,b,Aeq,beq,vlb,vub)

运算结果如下：

(2)例:投资的收益和风险

二、 基本假设和符号规定

一、问题提出

市场上有n 种资产i s （i =1,2，…，n ）可以选择，现用数额为M 的相当大的资金作一个时期的投资.这n 种资产在这一时期内购买i s 的平均收益率为i r ，风险损失率为i q ，投资越分散，

总的风险越小，总体风险可用投资的i s 中最大的一个风险来度量.

购买i s 时要付交易费，(费率

i p )，当购买额不超过给定值i u 时，交易费按购买i u 计算.

另外，假定同期银行存款利率是0r ，既无交易费又无风险.（0r =5%） 已知n =4

试给该公司设计一种投资组合方案，即用给定的资金M ，有选择地购买若干种资产或存银行生

息，使净收益尽可能大，使总体风险尽可能小.

基本假设：

1. 投资数额

M 相当大,为了便于计算，假设M =1； 2．投资越分散，总的风险越小； 3．总体风险用投资项目i s 中最大的一个风险来度量； 4．n 种资产i s 之间是相互独立的；

5．在投资的这一时期内, r i ,p i ,q i ，r 0为定值，不受意外因素影响； 6．净收益和总体风险只受 r i ,p i ,q i 影响，不受其他因素干扰.

符号规定： S i ——第i 种投资项目，如股票，债券

r i ,p i ,q i ----分别为S i 的平均收益率, 交易费率，风险损失率

u i ----S i 的交易定额 0r -------同期银行利率

x i -------投资项目S i 的资金 a -----投资风险度

Q ----总体收益 ΔQ ----总体收益的增量

三、模型的建立与分析

1.总体风险用所投资的Si 中最大的一个风险来衡量,即max{ qixi|i=1，2,…,n}

4. 模型简化：

四、模型1的求解

由于a 是任意给定的风险度，到底怎样给定没有一个准则，不同的投资者有不同的风险度.我们从a =0开始，以步长△a =0.001进行循环搜索，编制程序如下：

2．购买S i 所付交易费是一个分段函数,即

p i x i x i >u i 交易费 = p i u i x i ≤u i

而题目所给定的定值u i (单位:元)相对总投资M 很小, p i u i 更小,

可以忽略不计,这样购买S i 的净收益为(r i -p i )x i

3．要使净收益尽可能大，总体风险尽可能小,这是一个多目标规划模型: 目标函数 max

∑

=-n

i i i i x p r 0)( minmax{ q i x i }

约束条件 0(1)n

i i i p x =+∑

=M

x i ≥0 i =0,1,…,n

a ． 在实际投资中，投资者承受风险的程度不一样，若给定风险一个界限a ，使最大的一个

风险q i x i /M ≤a ，可找到相应的投资方案. 这样把多目标规划变成一个目标的线性规划. 模型1 固定风险水平，优化收益

目标函数： Q =max ∑

+=-11

)(n i i i i x p r

约束条件： M

x q i

i ≤a ∑

=+M x p i i )1(， x i ≥ 0 i =0,1,…, n

b ．若投资者希望总盈利至少达到水平k 以上，在风险最小的情况下寻找相应的投资组合.

模型2 固定盈利水平，极小化风险

目标函数： R = min{max{ q i x i }} 约束条件：∑

=-n i i i i x p r 0

)(≥k ， ∑=+M x p i i )1( ， x i ≥ 0 i =0，1，…，n

模型1为: minf = (-0.05, -0.27, -0.19, -0.185, -0.185) (x 0 x 1 x 2 x 3 x 4 ) T

x 0 + 1.01x 1 + 1.02x 2 +1.045x 3 +1.065x 4 =1 0.025x 1 ≤a

0.015x 2 ≤a s.t. 0.055x 3 ≤a

0.026x 4≤a x i ≥0 (i = 0，1， (4)

a=0;

while(1.1-a)>1

c=[-0.05 -0.27 -0.19 -0.185 -0.185];

Aeq=[1 1.01 1.02 1.045 1.065]; beq=[1];

A=[0 0.025 0 0 0;0 0 0.015 0 0;0 0 0 0.055 0;0 0 0 0 0.026];

b=[a;a;a;a];

vlb=[0,0,0,0,0];vub=[];

[x,val]=linprog(c,A,b,Aeq,beq,vlb,vub);

a

x=x'

Q=-val

plot(a,Q,'.')，axis([0 0.1 0 0.5])，hold on

a=a+0.001;

end

xlabel('a'),ylabel('Q')

结果：

a = 0.0030 x = 0.4949 0.1200 0.2000 0.0545 0.1154 Q = 0.1266 a = 0.0060 x = 0 0.2400 0.4000 0.1091 0.2212 Q = 0.2019 a = 0.0080 x = 0.0000 0.3200 0.5333 0.1271 0.0000 Q = 0.2112 a = 0.0100 x = 0 0.4000 0.5843 0 0 Q =0.2190

a = 0.0200 x = 0 0.8000 0.1882 0 0 Q =0.2518

a = 0.0400 x = 0.0000 0.9901 0.0000 0 0 Q =0.2673

四、实验三：插值

1.实验目的：

了解插值的基本内容．

2.实验内容：

例：

用分段线性插值法求插值,并观察插值误差.

（1）在[-6,6]中平均选取5个点作插值(xch11)

（2）在[-6,6]中平均选取11个点作插值(xch12)

6

6

,

1

1

)

(

2

≤

≤

-

+

=x

x

x

g

（3）在[-6,6]中平均选取21个点作插值(xch13)

（4）在[-6,6]中平均选取41个点作插值(xch14)

五、实验四：拟合

1.实验目的：

直观了解拟合基本内容．

掌握用数学软件求解拟合问题.

2.实验内容：

例: 对下面一组数据作二次多项式拟合

输入指令： x=0:0.1:1;

y=[-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2]; A=polyfit(x,y,2) z=polyval(A,x); plot(x,y,'k+',x,z,'r')

结果：

0317

.01293.208108.9)(2-+-=x x x f

六、MATLAB 在材料力学里的应用

1、实验目的：

通过构建可视化窗口，实现对平面问题：材料数据的直观、简便处理，得到需要的结论，具体功能如下：1.已知应变场，输入窗口后，能够得到应力场；2.已知应力场，根据材料的物性方程，判断材料是否发生屈服；3.已知某两点应力场，通过矩阵运算得到应力张量不变量，从而判断这两点是否处于同一应力状态。

数学背景与材料力学解释：

1.上述各类场量都易以矩阵形式表达，并且这些场量相互之间存在以下关系：

Mises 屈服准则：

平面应力物理方程：

斜方向应变公式：

主应力公式：

2.通过单向拉伸与扭转实验，可以测得弹性模量（E ）、剪切模量（G ）、泊松

)(1

y x x E

μσσε-=)

(1

x y y E

μσσε-=

G

xy

xy τγ=

?????++=++=++=xy y x N xy y x N xy y x N m l m l m l m l m l m l γ

εεεγεεεγεεε3323233

222

22221121211

比（u ）、应变场量，通过这些实验数据进行运用求解，是具有实际意义的。

求解过程：

在单向拉伸实验样品上按一定角度（Lij ）贴应变片（默认0度、45度、90度）,得到弹性模量（E ）、应变量（1N ε、2N ε、3N ε）。通过斜方向应变公式得到应变张量，再由物理方程得到应力张量。利用求出的应力场由主应力公式求得主应力，再由等效应力公式、Mises 屈服准则得到等效应力以及强度校核。

2、实验内容：

（1）.功能一：“开始运算”求应力应变的实现

function pushbutton1_Callback(hObject, eventdata, handles)

% hObject handle to pushbutton1 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) e1=str2num(get(handles.edit10,'String'))

e2=str2num(get(handles.edit11,'String'))

e3=str2num(get(handles.edit12,'String'))

u=str2num(get(handles.edit20,'String'))

E=str2num(get(handles.edit1,'String'))

res=2*e3-e1-e2

set(handles.edit16,'String',res)

set(handles.edit14,'String',e1)

set(handles.edit15,'String',e2)

G=E/(2*(1+u))

x=(u*E*e2+E*e1)/(1-u.^2)

y=(x-E*e1)/u

z=G*res

set(handles.edit17,'String',x)

set(handles.edit18,'String',y)

set(handles.edit19,'String',z)

x1=(x+y)/2+sqrt(((x-y)/2).^2+z)

y1=(x+y)/2-sqrt(((x-y)/2).^2+z)

z1=x+y+z-x1-y1

set(handles.edit21,'String',x1)

set(handles.edit26,'String',y1)

set(handles.edit27,'String',z1)

eq=sqrt(((x1-y1).^2+(y1-z1).^2+(z1-x1).^2)/2)

set(handles.edit29,'String',eq)

（2）.功能二：“强度校核”功能的实现

function pushbutton3_Callback(hObject, eventdata, handles)

% hObject handle to pushbutton3 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) s=str2num(get(handles.edit31,'String'))

e=str2num(get(handles.edit29,'String'))

if e

set(handles.edit30,'String','材料未发生屈服')

else

set(handles.edit30,'String','材料发生屈服')

end

功能三：“Clear ALL”的实现

% --- Executes on button press in pushbutton2.

function pushbutton2_Callback(hObject, eventdata, handles)

% hObject handle to pushbutton2 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) set(handles.edit10,'String','')

set(handles.edit11,'String','')

set(handles.edit12,'String','')

set(handles.edit20,'String','')

set(handles.edit1,'String','')

set(handles.edit16,'String','')

set(handles.edit14,'String','')

set(handles.edit15,'String','')

set(handles.edit17,'String','')

set(handles.edit18,'String','')

set(handles.edit19,'String','')

set(handles.edit21,'String','')

set(handles.edit26,'String','')

set(handles.edit27,'String','')

set(handles.edit29,'String','')

set(handles.edit30,'String','')

set(handles.edit31,'String','')

3.总结与评价：

通过利用材料力学中的斜方向应变公式、物理公式、主应力公式、Mises屈服等数学概念，运用MATLAB中的GUI可视化窗口模型构建，实现了直观、简洁的处理材料力学相关实验数据的目的，输入基础实验数据，便得到了有利于材料结构设计所需要的基本力学参量、也能对材料力学性能的评价起到重要作用。十分具有实际意义。

问题在于，这里处理的是平面问题，尤其是平面应力问题，可以推广到三维应力应变问题；而且斜方向应变公式也是取用的0度、45度、90度特殊方向。对于学习、验证某些普通应变情况不适用，但对于实验，我们可以认为控制这些变量，如：采用一定长宽比的长薄板，控制为平面应力状态，且由圣维南原理知薄板中部均匀受力；按特殊方向帖应变片。

七、MATLAB创建2048小游戏

1游戏规则：

1）.通过按键（W，A，S，D或方向键）控制小块的移动方向，当有两个小块相同时，把两者合并，产生一个新块（其值90%为2，10%为4），给予一定的分数，并将分数记录在数据文件中，显示在界面上。

2）.当最大块的值达到2048时，您将获得胜利，胜利后您可以选择继续游戏（当最大块的值为4096时，您将获得连胜，依次类推），或者退出游戏（下次游戏时从退出时的界面开始）。

3）.但是，如果小块不能再移动，您将输掉该游戏，此时，您可以选择新游戏或者退出游戏（下次游戏时重新开始）。

4）.在游戏过程中，您都有权停止游戏（按Esc键或者直接关闭界面），下次打开该游戏时将回到本次停止的状态，以便您继续游戏，当然你也可以随时新建游戏。5）.该游戏为了增加游戏的观赏性，在工具栏内特为您增加了几项设置，如新建游戏、游戏界面截图、动画、声音、演示、游戏得分记录等，您可以按照您的喜好进行设置。

6）.游戏还增加了数据分析的功能，您可以查看当前统计，也可以查看历史记录（只要您没有删除他们）。

7）.为了降低游戏难度，本游戏还提供了提示功能，按“空格”键即可提示，但连续按下空格键的时间间隔不得少于1s，否则会卡住，影响娱乐效果。

8）.游戏还提供了演示功能，您可以设置电脑演示步数，让电脑替您玩游戏。大概960步可以完成2048(只要小块可以一直移动)。

9）.由于AI提示程序在个别情况下无法得出最佳的移动策略，在这些情况下，您可以选择工具栏中“炸弹”按钮，用鼠标左键选中一个块将其炸掉（您无法在游戏开始时炸掉界面内唯一的块），然后继续您的游戏，当您想取消炸弹时，可以在游戏界面中点击鼠标右键。作为炸掉有数字的块的惩罚，您将损失该块的分数。原则上您可以炸毁界面内所有的数字块，但若界面内的块均被您炸掉，则游戏终止（除非您再次新建游戏）。

10）.点击“炸弹”按钮后，您也可以继续您的游戏，移动几步后再进行炸毁小块，比

如AI演示时进行炸毁操作。如果该小块没有数值，则炸毁无效，不扣分。若您出现无法移动的情况，在您使用“炸弹”按钮前，请不要在弹出的对话框中选择“新游戏”或“退出”，直接点击对话框红色“×”号关闭对话框后再使用炸毁功能，否则游戏将结束或关闭，而无法继续该局游戏。

2、游戏代码与运行结果

（1）.构建主坐标轴

handles.hAx = axes(...

'Parent', handles.hFig, ... % 父对象为handles.hFig

'Units', 'Pixels', ... % 单位像素

'Position', [10 10 400 400], ... % 位置

'Color', [187 173 160]/255, ... % 背景颜色

'XLim', [0 4], ... % X轴范围

'YLim', [0 4], ... % Y轴范围

'XTick', [], ... % X轴刻度为空

'YTick', [], ... % Y轴刻度为空

'XColor', [187 173 160]/255, ... % X轴颜色

'YColor', [187 173 160]/255, ... % Y轴颜色

'Box', 'on', ... % 显示边框

'XTickLabel', [], ... % X轴刻度标签为空

'YTickLabel', [], ... % Y轴刻度标签为空

'PlotBoxAspectRatio', [1 1 1]);

handles.hScore = uicontrol(...

'Parent', handles.hFig, ... % 父对象为handles.hFig

'Units', 'Pixels', ... % 单位像素

'Position', [10 415 200 50], ... % 位置

'Style', 'Text', ... % 类型为text

'FontUnits', 'Pixels', ... % 字体单位像素

'FontSize', 14, ... % 字体大小14号

'HorizontalAlignment', 'left', ... % 对齐方式左对齐

'BackgroundColor', get(handles.hFig, 'Color'), ... % 背景颜色

'String', sprintf('Score: 0\nBest: %d', handles.Best)); % 显示文本

handles.hMoves = uicontrol(...

'Parent', handles.hFig, ... % 父对象为handles.hFig

'Units', 'Pixels', ... % 单位为像素

'Position', [300 415 120 50], ... % 位置

'Style', 'Text', ... % 类型为text

'FontUnits', 'Pixels', ... % 字体单位像素

'FontSize', 14, ... % 字体大小14号