MATLAB(01)

今天开始写读书笔记，第一本书是《高等应用数学问题的MATLAB求解》。

MATLAB语言是科学研究者首选的计算机数学语言。大多数学生缺乏对应用数学问题的全面了解，不清楚什么问题能用数学描述，什么样的数学问题能求解，以至于走了很多弯路。通过学习MATLAB语言可以使数学问题求解能力显著提高，即使对某些数学公式理解不够深刻，只要学习了MATLAB语言也能容易地求解类似问题，毕竟科学研究者和工程技术人员通常不考虑为什么这样做，而只考虑结果是什么。借助MATLAB数学基础不深厚的人同样能轻易利用计算机解决高深的应用数学问题。

第1章 计算机数学语言概述

1.1为什么学习计算机数学语言

并不是所有数学问题都能手工推导的，手工推导有时是不易、不可靠、不精确甚至是不可能的，所以需要计算机的帮助。用计算机的方式有两种，其一是数值分析，其二是计算机语言，前者只能解数值计算问题，后者还可以解决像公式推导这样的问题。

许多专门的学科在介绍原理与方法时一般采用简单的例子，回避高阶和复杂的例子，这样得到的方法和结论只是理论上的，用于解决现实的实际问题往往是比较困难的，所以需要借助计算机数学语言。

1.2 数学问题的解析解与数值解

数学家往往关注解析解以及解的存在性的严格证明，而工程技术人员通常只关心最后的数值解。数学上，解析解不存在的情况很普遍，例如圆周率就没有解析解，工程技术人员也不关心圆周率是怎么求得的，只是想知道能满足精度要求的数值，一般的计算取3.1416即可，粗略的计算取3.14也未尝不可。

第2章 MATLAB语言程序设计基础

MATLAB的主要特点：

简介高效 一条语句顶C/C++数百条

科学运算能力强大 MATLAB以矩阵为基本单元

绘图功能 论文中必须给出二维或三维图形

工具箱 为各专业量身定制

动态系统仿真 转用ANSYS等专业软件

2.1 MATLAB程序设计语言基础

2.1.1 MATLAB语言的变量和常量

变量必须由字母打头，且区分大小写。

MATLAB保留的变量：

eps——浮点运算误差限，默认2.2201e-16，若某个值小于eps则会被视为0

i,j——循环变量，用于纯虚数量，应恢复i=sqrt(-1)

inf——无穷大，-inf为无穷小，若除数为0，MATLAB也不会中止运算，而是给出一个“除0”警告，并将结果赋成inf

NaN——not a number不定式，用于0/0，inf/inf运算

lasterr、lastwarn、pi

2.1.2 数据结构

数值型数据：一般为double()，在图形处理等计算中通常使用uint8()uint16()int8()int32()等，u表示无正负号，数字表示位数。

符号型数据：用于推导公式和解析解法，使用前首先需要将采用的变量申明为符号变量

syms list props，例如syms a b，props可以是real,positive等。

变精度算法函数vpa()可以将符号型数值以任意精度显示出来，vpa(pi)将显示圆周率的32位有效数字，

vpa(pi,100)将显示圆周率的100位有效数字，pi也可以换成矩阵。

其他数据结构：字符串型数据、多位数组、单元数组、类于对象。

2.1.3 MATLAB的基本语句结构

分两种：直接赋值语句、函数调用语句

直接赋值语句：

末尾带分号则不会显示运算结果，不带则直接显示

如果没有指定变量，则结果赋予保留变量ans

注释以%打头

矩阵输入：A=[1,2,3;4 5,6;7 8 9]

函数调用语句:[返回变量列表]=函数名(输入变量列表)

一个函数名对应一个.m文件或MATLAB内部函数。

[]可以用逗号和空格分隔[U,S,V]，(X,Y,Z)只能用逗号分隔，为了避免混乱以后所有分隔用逗号完成。

2.1.4 冒号表达式与子矩阵提取

v1=1:0.2:pi

向量v1=1 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0

v2=1:pi

v2=1 2 3

v3=pi:-1:1

v3=3.1416 2.1416 1.1416

v4=1:-1:pi

Empty matrix

提取子矩阵：

:表示整行或整列

end表示最后一行或一列

B1=A(1:2:end,[2,3,4]) %奇数行234列

B2=A([1,2,3],[2,3,4]) 3行234列

B1=A(end:-1:1,:) %上下翻转，最后一行排在第一行，列排序不变

2.2 基本数学运算

2.2.1 矩阵的代数运算

矩阵转置：如果矩阵元素含有复数元素，则转置后取其共轭复数值，称为Hermit转置，记为B=A*,MATLAB中用A’表示。而矩阵的一般转置用A.’表示。

加减法运算：A和B的维数若相同，则相应元素相加减；若二者之一为标量，则遍加减于每个元素；否则报错。

矩阵乘法：C=A*B，并不需要指定二者的维数，若不相容则报错

矩阵的左除：A\B表示A-1B

矩阵的右除：B/A表示B A-1

矩阵翻转：fliplr(A)左右翻转，flipud(A)上下翻转，rot90(A)逆时针翻转90度

矩阵乘方运算：只有方阵可进行乘方运算，B=A^x

点运算：这是一种特殊且非常重要的运算，即对应元素的直接运算，这要求二者维数相同。A.*B、A.^x、A.^A

2.2.2 矩阵的逻辑运算

逻辑运算的结果是一个矩阵，矩阵的元素非0即1。

与运算：C=A&B，只有当AB所有元素非0，C才为1，否则C取0

或运算：C=A|B，只有当AB所有元素为0，C才为0，否则C取1

非运算：C=~A，只有当A元素为0，C为1，否则C取0

异或运算：C=xor(A,B)，只有当AB元素一个非0，一个为0，C才为1，否则取0，即只有二者逻辑结果不同才可取1，若相同则取0

2.2.3 矩阵的比较运算

一般比较运算的结果是一个矩阵，矩阵的元素非0即1。

C=A>B仅当A的元素大于B的对应元素C才取1，否则C取0.

类似的还有< 、>=、<

=、==、~=等

find()、all()、any()函数

find()函数用于从矩阵中提取满足条件的元素的下标，下标的表示方式有两种：

第一种，依次将矩阵的第一列、第二列•••排列成一个列向量，输出结果为这个列向量的下标，例如A=[1,2,3;4,5,6;7,8,9]变形为列向量就是A1={1,4,7,2,5,8,3,7,9}，find(A>=5)’的结果是3

5 6 8 9即A1列向量中的这些下标的元素满足大于等于5的筛选条件。

第二种，更为直观地输出行列下标，[i,j]=find(B>4);[i,j]的输出结果为

ans =

3 1

2 2

3 2

2 3

3 3

左侧为行标，右侧为列标

all()、any()函数返回的都是一个行向量，元素的数量跟A矩阵的列数相等，每个元素非0即1。

all(A>4)返回的行向量表示A矩阵中哪一列的所有元素都大于4

any(A>4)返回的行向量表示A矩阵中哪一列中存在至少一个元素大于4

2.2.4 解析结果的化简和变换

变量s用于表示解析式之前必须首先声明符号变量syms s

s1=simple(s)自动选择最简格式进行化简

[s1,var1]=simple(s)，字符串变量var1输出自动选择的化简方法的名称

常用的化简方法有：

collect()合并同类项

expand()展开多项式

factor()因式分解

numden()提取多项式的分子和分母

sincos()三角函数化简

subs()用于变量替换

>> syms a b c x m p q

>> f=a*x^2+b*x+c;

>> g=subs(f,a,m)

g =m*x^2 + b*x + c

>> h=subs(g,{b,c},{p,q})

h =m*x^2 + p*x + q

2.2.5 基本数论运算

floor()向负无穷方向取整

ceil()向正无穷方向取整

round就近取整

fix趋零取整

[n,m]=rat(A)将A中所有元素变换成最简有理数，分子存储在矩阵n中，分母m中

rem(A,B)A中元素对B中元素求模得出的余数

gcd(n,m) 最大公约数

lcm(n,m)最小公倍数

factor质因数分解

isprime是否为质数，以0/1表示

A=1:100;B=A(isprime(A))

2.3 流程结构

;表示当前代码未结束输入，但是需要换行，所以并不会立刻执行计算；有时一行代码中间也使用;表示在一行中书写原本应该多行书写的代码。

for、while、if、switch case、try catch

如果能对整个矩阵进行运算时，尽量不要采用循环结构

s=0;

for i=1:100

s=s+1;

…

end

case{,,,}

try,A

catch,B

end

首先试探性地执行A，如果没有错误则end，否则将错误信息赋值给lasterr并中止A转而执行B。A通常是高效但不稳定的算法，B通常是低效但稳定的算法，B也可用于说明A中的错误原因

2.4 函数编写与调试

2.4.1 基本结构

两种源程序文件格式：M-文件、M-函数。常用后者。

保留变量nargin和nargout用于确定函数的输入和输出变量的个数。

[输出变量列表]=函数名(输入变量列表)，均用逗号分隔。

help 函数名：显示注释语句内容，但空格行下方的注释不被显示。

if abs(n-floor(n))>eps|n<0

error(‘n应该是非负整数’);

上面的代码涉及abs、floor、eps、error、|等符号，用于判定是否为非负整数。

函数可以递归调用：

function srhf_fact(n)

if n>1

k=n*srhf_fact(n-1);

elseif any[0,1]==n

k=1

end

2.4.2 可变输入输出个数

如何建立起无限个输入或输出变量的函数？这种方式被广泛采用。

conv(A,B)用于求多项式的乘积，多个多项式的乘积只能嵌套使用此函数，即conv(C,

conv(A,B))。这样很不方便，可以考虑编写一个新的函数convs，其输入变量不限于2个，而是无限个，在MATLAB中“无限”并不是纯粹数学上的无限，而是指输入或输出变量的个数不再于函数编写时加以限制，所以在函数中需要知道用户在输入变量列表中给出了几个变量，这需要调用length函数。

函数输入变量列表可以是一个变量(实际上是一个向量)

varargin，此变量是保留变量，它表示多个变量的集合，那么length(varargin)的输出变量是一个正整数，表示varargin中具体包含的变量数量。

function a=convs(varargin)

a=1;

varnum= length(varargin);

for i=1:varnum

a=conv(a,varargin{i});

end

end

2.4.3 inline函数与匿名函数

一个简单的函数，没有必要单独创建一个.m文件，只用一行代码进行描述即可。

function a=f(x,y)

a=sin(x.^2+y.^2)

end

上面的函数是在.m文件中的代码，而有些情况下，如微分方程求解和最优化问题，这样的函数可以采用inline函数的形式给出。

f=inline(‘sin(x.^2+y.^2)’,’x’,’y’)

以上两种函数定义形式，无论采用哪种，都可以通过z=f(3,6)直接求出z=0.8509，显然inline函数更为简便。

匿名函数跟inline函数类似，区别有二：更简洁、可直接使用工作空间中的变量。

简洁体现在f=@(x,y) sin(x.^2+y.^2)

变量的优势十分明显：一般函数中的变量，无论是输入变量还是输出变量只在函数内部使用，跟工作空间中的变量是相对隔离的，工作变量作为输入变量代入函数，工作变量便替换函数变量进行计算，最后函数将输出变量的结果代入工作变量返回工作空间。这样做虽然严谨，有效避免了变量使用上的混乱，但有时不够方便，匿名函数的引入，有效地将工作变量和函数变量结合起来，使用上更为灵活，但有一点需要注意，匿名函数直接引用工作变量，只是引用当时的工作变量的值，引用完成之后，如果工作变量的值发生变化，匿名函数所引用的工作变量值不会随着变化，仍然是之前的值。

2.5 二维图形绘制

2.5.1基本函数：plot(t,y)

t1,t2,t3时刻的函数值分别是y(t1),y(t2),y(t3)，将其输入MATLAB就是两个向量t=[t1,t2,t3]，y=[

y(t1),y(t2),y(t3)

]，这种曲线实际上是各个时刻对应函数值之间的折线，当t向量足够密集，看起来就是曲线了，所以在某些特殊区间(慢变化区间)可以取点较稀疏，以提高运算速度，在另一些区间(快变化区间)则需要较密集，以得到这些特殊区间的较为准确的图形。

(1) t、y均为向量，如上，还有更为简单的形式plot([t1,t2,t3],[y1,y2,y3])

(2)

t为向量，y为矩阵。将y矩阵的每一行作为y向量与t共同绘制一条曲线，曲线数量跟y矩阵行数相等，这要求y矩阵的列数跟t向量的元素个数相等。

(3) t、y均为矩阵。要求t、y矩阵行列数均相等，取t矩阵的一行与y矩阵的同一行组成向量对进行绘制。

(4)

以上三种组合形式，如果不止一对t、y，而是多对，那么不必多次使用plot函数，可以这样一次绘出所有的图形plot(t1,y1,t2,y2,t3,y3)或plot(t1,y1,Option1,t2,y2,

Option2,t3,y3, Option3)，Option的具体用法参见P30.

图形的附属命令：

grid on/off：打开或关闭网格

hold on/off：保护当前坐标，以后再用plot函数时新的曲线将叠印在现有坐标系上

title()：添加标题

xlabel()和ylabel()：给坐标轴加标注

(abs(x)>n)：将此式与其他解析式点乘，当x在(-n,n)区间以外此式取1，否则取0，这是分段函数在MATLAB中的表达方式。

所得图形中有多个对象，如曲线、坐标轴、图形窗口等，每个对象都可以设置具体属性，属性也可以被读取。

set(句柄,’属性名1’,属性值1,’属性名2’,属性值2…)

v=get(句柄,’属性名1’)

属性名必须加单引号，属性值不加。

2.5.2 特殊意义的图形绘制函数：参见P32.

2.5.3隐函数的绘制

隐函数：f(x,y)=0

ezplot(‘f(x,y)’)或ezplot(‘f(x,y)’,[-m,n]) 仅在[-m,n]上绘制隐函数图形

2.6 三维图形绘制

2.6.1三维曲线绘制

plot3(x,y,z)或plot3(x1,y1,z1,选项1,…)

2.6.2 三维曲面绘制

已知z=f(x,y)，先用meshgrid()绘制网格矩阵数据x,y，再基于此平面绘制三位曲面，mesh()网格图，surf()表面图，surfc()带等高线，surfl()带光照，contour()等高线，contour3()三维等高线。例：

[x,y]=meshgrid(0:10,-5:5) %在(x,y)平面上绘制一个正方形网格矩阵

z=x.^2+y.^2

surf(x,y,z)

axis(-31,31,0,62,0,1) %将所得图形的立方体内的部分放大，名为重设坐标系，增大可读性

shading faceted/interp/flat %修饰显示形式，默认第一个

set(gca,’xlim,[-5,5]) %将x轴限制在一个区间内

2.6.3 三位图形的视角

view(a,b)，方位角a表示从-y轴向+x轴方向逆时向旋转的角度，仰角b表示从xy平面向+z轴方向旋转的角度，默认情况下a=-37.5，b=30.

下图是我用MATLAB绘制出的第一张图，纪念一下。

z=x.^2+y.^2