【数学建模】Matlab 编程

        MATLAB是美国MathWorks公司自20世纪80年代中期推出的数学软件，具有优秀的数值计算能力和卓越的数据可视化。由于Maltab编程方便，有大量内部函数和工具箱可以使用，作图也
十分方便，因此在数学实验和数学建模竞赛中，我们常使用Matlab作为我们的编程工具。

1. 常用函数介绍

(1) 三角函数
sin-----正弦           sinh-----双曲正弦        asin-----反正弦        asinh----反双曲正弦
cos-----余弦          cosh-----双曲余弦        acos-----反余弦        acosh----反双曲余弦
tan-----正切          tanh-----双曲正切        atan-----反正切        atanh----反双曲正切

(2) 指数函数与对数函数
exp-------指数        log--------e为底的对数         log10-----常用对数      sqrt-----平方根
与复数有关的函数
abs------模或绝对值      angle-----幅角          conj-----复共轭
imag----虚部          real------实部

(3) 含入函数及其它数值函数
fix-------向0舍入      floor----向负无穷舍入        ceil-------向正无穷舍入
round---四舍五入      rem(a,b)---计算a/b的余数

如fix(2.3)=2   floor(2.3)=2   ceil(2.3)=3   round(2.3)=2

如fix(-2.3)=-2   floor(-2.3)=-3   ceil(-2.3)=-2   round(-2.3)=-2
 

(4) 有关向量的函数

min(x): 向量x的元素的最小值               max(x): 向量x的元素的最大值
mean(x): 向量x的元素的平均值            median(x): 向量x的元素的中位数
std(x): 向量x的元素的标准差                 diff(x): 向量x的相邻元素的差
sort(x): 对向量x的元素进行排序            length(x): 向量x的元素个数
norm(x): 向量x的Euclidean长度            sum(x): 向量x的元素总和
prod(x): 向量x的元素连乘积                  cumsum(x): 向量x的累计元素总和
dot(x, y): 向量x和y的内积                      cross(x, y): 向量x和y的外积
cumprod(x): 向量x的累计元素总乘积

2. 矩阵常见计算

(1) 矩阵输入
最简单的方法是把矩阵的元素直接排列在方括号中。每行内的元素间用空格或逗号隔开，行与行之间用分号隔开。
如：

A = [1 4 7
       3 6 9
       6 7 4]

输入为

A=[1,4,7;3,6,9;6,7,4]
A=[1,4,7;
   3,6,9;          
   6,7,4]

输出结果为: 
1 4 7
3 6 9
6 7 4

(2) 矩阵的转置

矩阵的转置用符号'来表示。

如：

A=[1,4,7;3,6,9;6,7,4];
B=A'

也可直接转置：

B=[1,4,7;3,6,9;6,7,4]'

则显示为：

B=
1 3 6
4 6 7
7 9 4

(3) 矩阵的加减
矩阵的加减使用的是“+”和“-”运算符。进行矩阵加减运算必须是同型矩阵。
如：

A=[1,3,6;  
   4,5,7;
   7,8,9];
B=[3,5,7;  
   2,4,6;
   1,3,9];
C=A+B     

显示结果为：
C = 4 8 13
       6 9 13
       8 11 18

(4) 矩阵乘法

矩阵乘法用符号*表示。要求前一矩阵列数与后一矩阵行数相同。
如：

A=[1,4,7;  
   2,5,8];
B=[4,5,9;  
   1,7,8;
   3,2,1];
C=A*B     

结果为：
C = 29 47 48
       37 61 66

(5) 矩阵与数的乘法
如：

A=[1,5,8;
   2,6,9];
B=3*A

结果为:
B = 3 15 24
       6 18 27

(5) 矩阵的行列式

求方阵A行列式,用det(A)表示
如：

A=[1.3,6;
   2,5,8;
   3,9,11];
Z=det(A)

则结果为：

Z=7

(6) 矩阵求逆

非奇异矩阵A求逆用inv(A)表示。
如：

A=[1,3,6;
   2,5,8;
   3,9,11];
Z=inv(A)

结果为:
Z =-2.4286   3.0000   -0.8571
      0.2857  -1.0000   0.5714
      0.4286   0              -0.1429

如果验证，可计算C=A*Z

C= 1.0000    0             -0.0000

     -0.0000   1.0000     -0.0000

     0              0              1.0000

利用逆矩阵可以解方程组。

如：AX=b

A=[1,3,6;2,5,8;3,9,11];
b=[3,6,7]'
X=inv(A)*b

结果为:
X = 4.7143
    -1.1429
     0.2857

或用X=A\b也可求解。
 

X=A\b还可以求解齐次方程组。

下面进行函数拟合：如因变量y与自变量x之间存在关系: y = a + b * e^(-x)
观测数据对为：
x    0.0   0.3   0.8   1.1   1.6   2.3
y   0.82  0.72  0.63  0.60  0.55  0.50

由此可建立齐次方程组  AX = Y 其中

即：

可建立Matlab的m文件

t=[0.0  0.3  0.8  1.1  1.6  2.3]'
Y=[0.82  0.72  0.63  0.60  0.55  0.5]'
A=[ones(size(t)),exp(-t)]
X=inv(A'*A)*A'*Y
%X=A\Y

可得结果为:
X = 0.4760  0.3413 即a=0.4760, b=0.3413

函数拟合为 y = 0.476 + 0.3413 * e^(-x)

Matlab的M文件如下

t=[0.0  0.3  0.8  1.1  1.6  2.3]'
Y=[0.82  0.72  0.63  0.60  0.55  0.5]'
A=[ones(size(t)),exp(-t)]
x=A\Y
n=500; %数据点数
tt=zeros(n,1);
yy=zeros(n,1);
dt=2.3/n;
for i=1:n
    tt(i)=i*dt;
    yy(i)=x(1)+x(2)*exp(-tt(i));
end
plot(t,y,'*b',tt,yy,'r')

(7) 矩阵特征值

如果A为方阵，满足AX = λX 的λ称为A的特征值，X称为A的特征向量。计算A的特征值用eig(A)表示。
如

A=[1,3,6;2,5,8;3,6,8];
Z=eig(A)

结果为:
Z = 15.2382
     -1.3365
      0.0982
 

3. 函数作图

(1) 二维平面曲线作图函数 plot(x,y,'s')
其中x和y是长度相同的向量，s表示线型和颜色。如果作多条曲线在同一图上，则用函数
plot(x1,y1,'s1',x2,y2,'s2',...,xn,yn,'sn')

如将 sin(x)和cos(x)同时作在一张图上，区间取[0,2π]，程序如下：

x=0:0.1:2*pi;
y1=sin(x);
y2=cos(x);
plot(x,y1,'r',x,y2,'b');

(2) 多窗口作图

将屏幕分为几个窗口分别作图，subplot(m,n,k) 表示将窗口分为 m × n 个，当前图在第k个窗口完成。如在第一个窗口作 y=sin(x)，第二个窗口作 y=cos(x)，第三个窗口作 y=√x， 第四个窗口作 y=ln(x)。程序如下：

x1=0:0.1:2*pi; y1=sin(x1);
x2=-pi:0.1:pi; y2=cos(x2);
x3=0:0.1:10; y3=sqrt(x3);
x4=2:0.1:10; y4=log(x4);

subplot(2,2,1); plot(x1,y1); title('y=sin(x)'); grid on
subplot(2,2,2); plot(x2,y2); title('y=cos(x)'); grid on
subplot(2,2,3); plot(x3,y3); title('y=sqrt(x)');grid on
subplot(2,2,4); plot(x4,y4); title('y=ln(x)'); grid on

(3) 直方图作图 hist

二维直方图，可以显示出数据的分布。

用法：count = hist(X) 把向量x中的元素放入等距的10个条形中，且返回每一个条形中的元素个数。

          count = hist(X,center) 参量X为向量,把X中元素放到 m (m=length(center)) 个由center中元素指定的位置为中心的直方图中。

          count = hist(X,number) 参量number为标量，用于指定条形的数目。

          [count,center] = hist(X) 返回向量X中包含频率计数的count与条形的位置向量center，可以用命令bar(center,count)画出条形直方图。

如作1000个服从正态N(10,25)数据的直方图，程序如下：

X=normrnd(10,5,1000,1);
hist(X);

采用下面的程序也可以：

X=normrnd(10,5,1000,1);
[count,center]=hist(X);
bar(center,count);

(4) 二维图形注释命令 grid

功能：给图形的坐标轴增加分隔线。该命令会对当前坐标轴的属性有影响。

用法：grid on 给当前的坐标轴增加分隔线。
           grid off 从当前的坐标轴中去掉分隔线。
           grid 转换分隔线的显示与否的状态。
           grid(axes_handle,on|off) 对指定的坐标轴

           axes_handle是否显示分隔线。

(5) 空间曲线作图

Matlab提供了三维曲线作图函数plot3。调用格式：plot3(x,y,z,'s')，其中x，y和z是长度相同的向量,s表示线型和颜色。作空间螺旋线：

t=0:0.01:8*pi;
x=cos(t);
y=sin(t);
z=t;
plot3(x,y,z,'r');

(6) 三维曲面作图

命令1 mesh
功能：生成由X，Y和Z指定的网线面
用法：mesh(X,Y,Z) 画出三维网格图

如作z=cosx.siny曲面图，程序如下：

[X,Y] = meshgrid(-3:0.1:3,-4:0.1:4);
Z=cos(X).*sin(Y);
mesh(X,Y,Z);
xlabel('x');
ylabel('y');
zlabel('z');

 4. 基本语句

(1) for语句

   for x=a:d:b
       (command)
   end

a为起始点，b为终止点，d为区间间隔。循环体内为执行语句。

求1+3+5+.....99。for语句程序 s=0;

for i=1:2:99
    s=s+i;
end
s

结果为2500

(2) while 语句

while expression
    (command)
end

用while语句实现如下：

s=0; i=1;
while i<100
    s=s+i;   i=i+2;
end
s

(3) if-else-end语句

该语句通常有三种形式
单分支

if expression
    (command) 
end     

双分支

if expression
    (command1)  
else
    (command2)   
end

多分支

if expression
    (command1)
elseif expression2
    (command2)
...
else
    (commandn)
end

如用该语句求1到100中被3整除余1的数和为s1，被3整除余2的数和为s2，被3整除的数和为s3。程序如下：

s1=0;s2=0; s3=0;
for i=1:100
   if(mod(i,3)==1) s1=s1+i;
   elseif (mod(i,3)==2) s2=s2+i;
   else s3=s3+i;  end
end
fprintf('s1=%3d s2=%3d s3=%3d\n',s1,s2,s3);

结果如下:
s1=1717  s2=1650  s3=1683

5. M文件

Matlab中M文件有两种。脚本M文件和函数M文件。

建立M文件的方法是：在Matlab的窗口中单击File➔New➔M-File。打开M-File窗口，在该窗口中输入程序。再以.m为扩展名存储。要运行M文件，F5或在命令窗口中输入该文件名，按确定键就可以了。

当我们要编写一个函数便于主程序调用时，该函数通过写为M文件。这种M文件称为函数M文件。函数M文件是文件名后缀为m的文件，这类文件的第一个非注释行必须以一个特殊字符function开始，格式为：function 函变量名=函数名(自变量名)

函数M文件和脚本M文件主要有以下差异：

(1) 函数M文件的文件名必须与函数名相同，脚本M文件则可以任意取合法的文件名。

(2) 脚本M文件没有输入参数与输出参数，而函数M文件有输入与输出参数，对函数进行调用时，可以按少于函数M文件规定的输入与输出变量个数，但不能多于函数M文件规定的输入与输出变量个数。

(3) 脚本M文件运行产生的所有变量都是全局变量，而函数M文件的所有变量除特别声明外都是局部变量。

如计算函数 z = (x^2 - 2x) * e^(-x^2 - y^2 - xy) 在(0.1,0.2)处的函数值。

编写M文件fun.m

function z=fun(x)
z=(x(1)^2-x(1)*2)*exp(-(x(1)^2-x(2)^2-x(1)*x(2)));
return;

在命令窗口中输入
>>x=[0.1,0.2]
>>z=fun(x)

得到结果为:
z=-0.1772