lecture+for+matlab

1 Matlab 在基础物理学中的应用 （内部讲义） 贾砚宾 管 靖 编 北京师范大学物理系 2004．10 2 Matlab 在基础物理学中的应用 这是选读内容，由教师组织学习当然更好，同学自学也是可以 的. 老师不一定是专家，韩愈说过：“师不必贤于弟子，弟子未必不 如师. ”“玩”计算机，正是年轻人的特长. 在这里，你有在大学一年 级就超越老师的机会. 计算机是“玩”会的，你拿出玩电子游戏的热情与专注，很快就 会成功. Matlab 在科研和 工程 路基工程安全技术交底工程项目施工成本控制工程量增项单年度零星工程技术标正投影法基本原理 设计中被广泛使用，它几乎可以做你 想做的所有事，学会了你可以终生受益，这比玩电子游戏有价值得 多. 现在，我们就可以借助它学习物理学. 你不要学了再做，要在做中去学. 装上软件，边看书边实践. 先 按书输入程序，让它能运行. 再模仿它，做一点点修改……发挥你 的创造性，做你想做的事…… 这里只是入门，入门后有很多书可看. 遇到问题，最快捷的方 法是向别人请教. Matlab 有很好的帮助系统（§7 之 3），但它是英文界面，开始 看不懂是有情可原的. Matlab 是 “Matrix Laboratory”的缩写，是一种功能强、 效率高、便于进行科学和工程计算的交互式软件包. 由于使用 Matlab 编程运算与人进行科学计算的思路 和 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 达方式完全一致，所以不象学习其它高级语言那样难 于掌握. 实践证明，学习者可在短短几个小时的学习和使 用中就能初步掌握 Matlab 的基础知识，从而使学习者能够 进行高效率和富有创造性的计算. Matlab 对于学习者的编程语言基础要求不高，库函数 和编程语句丰富多样且简单易学，在数据可视化上也有独 特的优势. 学习者不需要投入太多的时间在学习编程语言 知识上，可以直接利用软件提供的丰富的函数，编写较简 单的程序即可解决许多普通物理学的问题. 下面将结合实例，简要介绍 Matlab 的基础知识和在普 3 通物理学中的应用. 注意：带*的内容开始不必看，使用时再查阅. §1 Matlab 使用介绍 在 Windows 窗口中用鼠标双击 Matlab 图标 即可进 入 Matlab 的工作窗口（Command Window），如图 1 所示. 没 有图标可利用 Matlab\bin 目录下的 Matlab.exe 文件在桌面 上建立一个快捷方式. 退出 Matlab 的方法有三种：单击工作窗口右上角的关 闭按钮 ；用菜单 FileÆExit Matlab 命令；或者直接在工 作窗口中输入 quit 后回车. 工作窗口是标准的 Windows 窗口形式，用户在命令窗 口中输入各种指令，进行运算；在左侧的变量窗口中监控 当前所创立的所有变量. 当前工作路径 命令窗口 变量窗口 图 1 Matlab 工作窗口 Current Directory 是系统的当前工作路径，Matlab 对函 数或文件等进行搜索，用户每次文件的创建、保存都在这 个路径下进行. 初次启动 Matlab 时系统的默认工作路径是 Matlab目录下的Work子目录，如果要改变当前的工作路径， 4 可以单击如图 2 所示的路径栏右侧的 ，在弹出的路径选 择对话框内选择想要设置的路径. 图 2 工作路径栏和选择工作路径对话框 §2 Matlab 能做什么 本节先给出使用 Matlab 辅助普通物理学学习的几个简 单的例子，意在引导读者对 Matlab 的功能特点及语句编程 有一个概括的了解. （1）绘制简谐振动的振动曲线 简谐振动的运动学方程是 ( )ϕω += tAx sin ，可以根 据这个方程利用 Matlab 的绘图语句画出质点简谐振动的 位移曲线. 在命令窗口中输入：（lt21.m） a=1; omiga=1; %为变量赋值 t=0:0.05*pi:4*pi; %为变量赋值 phi=0.1*pi; %为变量赋值 x=a*sin(omiga*t+phi); %建立表达式并运算 plot(t,x) %以 t 为横坐标，x 为纵坐标绘图 运行结果如图 3 左图所示. 5 图 3 简谐运动的位移曲线 在这个例子中，a 是振幅 A，omiga 是圆频率ω，phi 是初相位ϕ，t 是自变量，开始为这些参量以及自变量赋 值. 可以看出，Matlab 不像其他编程语言那样必须进行变 量的预定义，创建变量和为变量赋值是同时完成的. 另外， 对时间 t 的赋值是一组数据，即从 0 到 4π每隔 0.05π 取一 个值，形成一个数组：t=[0 0.1571 0.3142 … 12.2522 12.4093 12.5664 ] （数组中的数用空格或逗号分开），总共有 81 个 元素. 语句中的 pi 是 Matlab 默认的常数 π. x=a*sin(omiga*t+phi)是建立表达式并运算，sin（）函 数是 Matlab 的内置函数，直接调用即可，由此得出对应每 一个 t 值的质点的位移 x，x 是一个与 t 所含元素个数相同 的数组. 最后的 plot(t,x)是二维绘图语句，表示以 t 为横坐标， x 为纵坐标绘图，得出的就是正弦曲线. 如果对各个参量取 不同的值，得出的曲线形状也就不一样. 图 3 右图是选择 不同参量时的曲线对比情况. 对于数组 t 的赋值，一般格式为 t=t0:tstep:tend. t0 为数 组中第一个元素的值，tend 为最后一个元素的值，tstep 为 取值的间隔. Matlab 会根据计算的结果自动调整坐标轴范围，当然 用户也可以通过语句对坐标轴进行控制. 6 （2）等量异号点电荷的电势分布 这个例子将介绍二维网格和三维曲面绘图的语句. 物 理情景是 Oxy 平面上在 x=2, y=0 处有一正电荷，x=-2, y=0 处有一负电荷，根据 r qU 04πε= 计算两点电荷电场中电势 的分布, 2020 )()( yyxxr −+−= . 在命令窗口中输入：（lt22.m） [x,y]=meshgrid(-5:0.2:5,-4:0.2:4); %建立数据网格 z=1./sqrt((x-2).^2+y.^2+0.01)-1./sqrt((x+2).^2+y.^2+0.01); %电势的表达式 mesh(x,y,z) %三维曲面绘图 运行结果如图 4 所示. 选定一系列的 x 和 y 后，就组成了平面上的网格点， 再计算对应每一点上的 z值. -5:0.2:5,-4:0.2:4分别是选取横 坐标与纵坐标的一系列数值，meshgrid 是生成数据网格的 命令， [x,y]是 xy 平面上的坐标网格点. z=1./sqrt((x-2).^2+y.^2+0.01)-1./sqrt((x+2).^2+y.^2+0.01) 是场点（x, y）的电势，其中 sqrt( )是 Matlab 默认的函数： 求变量的平方根. 当场点即在电荷处时，会出现分母为零 的情况，因此在 r 里加了一个小量 0.01，这样既可以完成 计算，又不会对结果的正确性造成太大影响. 图 4 等量异号点电荷的电势分布 图 5 λ=0.8 时的叠代情况 7 ) 另外需要注意的是表达式中的“./”、“.^”是对数组运算 的算符，含义与数值运算中的“/”、“^”相同，不同之处是 后者只对单个数值变量进行运算，而前者对整个数组变量 中的所有元素同时进行运算. mesh 是三维网格作图命令，mesh(x,y,z)画出了每一个 格点（x, y）上对应的 z 值（电势）. （3）Logistic 叠代图 ( nnn xxx −=+ 11 λ 是一个著名的非线性叠代方程，称 为Logistic叠代，给定一个初值x0（0k 则返回空值 j:i:k 等价于[j,j+i,j+2*i,……,k] i>0 则要求 jk，否则返回空值 此外，冒号还可以用作矩阵的下标，以及部分的选择 矩阵的元素，执行循环操作等，后面几节中会陆续介绍. ②百分号“％”此符号在命令行中表示注释，即在一行 中百分号后面的语句都被忽略而不被执行. ③连续点“…”如果一条命令很长，一行容不下，可以 用 3 个点加在一行的末尾，表示此行未完，而在下一行继 续. 例如在建立等量异号点电荷电场中电势的表达式时， 可以这样写： z=1./sqrt((x-2).^2+y.^2+0.01)… -1./sqrt((x+2).^2+y.^2+0.01); 这样虽然占了两行，但 Matlab 仍然认为这是一条语 句，不会出现错误. ④分号“；” 分号用在每行命令的结尾，要求执行命 令但不显示计算结果. *（3）关系运算符 关系运算符主要用于在数与数、矩阵与矩阵之间进行 比较，基本的有： 运算符 < > <= >= == ~= 功能 小于 大于 小于 等于 大于 等于 等于 不等 于 *（4）逻辑运算符 在 Matlab 中包含与“&”、或“|”、非“~”、异或“xor”4 种逻辑运算符. （三）变量与表达式 用运算符把数字、变量和函数组合在一起，就建立了 一个表达式. 例如前面例子中的 x0=abs(sin(randn))， x=a*sin(omiga*t+phi)等. 在 Matlab 中，一个变量可以通过 给它分配一个数值或表达式来定义，如下所示： variable =expression 一个变量的值可以通过输入它的名字=值（或 10 表达式），并按回车键获得，Matlab 以显示这个变量的名 字和值作为回答. 如果这个变量并不存在，就显示一个错 误信息，如图 6 所示. 图 6 显示变量 在 expreesion 之后可以加分号后按回车键，也可以直 接按回车键. 没有结尾分号的每个命令在屏幕上显示出其 结果；若结尾带分号，就执行计算，但计算结果并不显示. 如果不指定变量而直接输入 expression 项，则 Matlab 用 ans（ answer 的缩写）显示这个值. *（4）管理变量 当创建了很多不同类型的变量后，就需要获得当前工 作窗口中的变量的信息，并清除不再需要的变量. Matlab 提供了一些命令专门用来管理变量. 名称 含义 who,whos 列出当前工作窗口的变量 clear 从内存中清除变量和函数 size 获得矩阵的尺寸（矩阵的维数、元素个数等） length 获得矢量的长度 disp 在显示变量名称的同时显示变量的内容 *（5）基本数学函数 Matlab 中有许多已定义好的数学函数，用户可以直接 调用，例如前面例子中的 sin，randn，abs 等，常见的数学 函数可参考下表. 11 ①三角函数与双曲函数 名称 含义 名称 含义 名称 含义 sin 正弦 sec 正割 asinh 反双曲正弦 cos 余弦 csc 余割 acosh 反双曲余弦 tan 正切 asec 反正割 atanh 反双曲正切 cot 余切 acsc 反余割 acoth 反双曲余切 asin 反正弦 sinh 双曲正弦 sech 双曲正割 acos 反余弦 cosh 双曲余弦 csch 双曲余割 atan 反正切 tanh 双曲正切 asech 反双曲正割 acot 反余切 coth 双曲余切 acsch 反双曲余割 ②对数函数 名称 含义 名称 含义 名称 含义 Exp e 为底的指 数 log10 10 为底的对数 pow2 2 的幂 Log 自然对数 log2 2 为底的对数 sqrt 平方 根 ③复数函数 名称 含义 名称 含义 名称 含义 Abs 绝对值 conj 复数共轭 real 复数实部 Angle 相角 imag 复数虚部 ④统计函数 名称 含义 名称 含义 min 最小值 max 最大值 mean 平均值 median 中位数 sum 总和 diff 相邻元素的差 sort 排序 length 个数 12 §4 矩阵与数组 （一）定义矩阵 Matlab 中数据的基本格式是矩阵. 2 维矩阵是一个以 行和列排列的元素矩形表. 如果有 m 行、n 列，这个矩阵 的大小就是 m×n. 如一个 2×3 的矩阵： . ⎥⎦ ⎤⎢⎣ ⎡ 654 321 第 1 行是[1 2 3]，第 2 列是 . ⎥⎦ ⎤⎢⎣ ⎡ 5 2 这里之所以说到矩阵，是因为 Matlab 就是矩阵试验室. 实际 上读者没学过矩阵也没关系，等你学过以后，Matlab 会显示出更强 大的功能. 现在不必去学矩阵，我们先把矩阵理解为或一行、或一 列、或 m 行 n 列的数组，我们也主要只进行数组的运算. 矢量在三 维空间里有三个分量，现在我们把矢量推广为有多个分量；如果数 组只有一行，它就是一个行矢量；如果数组只有一列，它就是一个 列矢量. 下面我们就使用矩阵、矢量进行讲解，你把它想成数组就 行了. 矩阵的元素，即数aij，通常是实数，但也可以是复数. aij是指第i行、第j列的数. 在上例中，有a21= 4. 当矩阵仅由一行组成时，它是一个特例，就是一个行 矢量. 如果矩阵仅有一列，就是一个列矢量. 矢量中元素的 数量就是矢量的长度. 如果矩阵的维数是 1，它是一个标 量，就是一个数. 二维矩阵的建立可以有多种方法实现，最简单的方法 是由方括号[ ]包围的逐行给定元素. 相同行中的元素用空 格‘’或逗号‘,’分隔，不同行由分号‘;’或回车键分隔，如图 7 所示. 如果定义一个标量，则方括号就不需要了. 13 图 7 矩阵建立和元素的引用 （二）矩阵元素的标识与修改 一个矩阵 A（或一个矢量 A）的指定元素用 A(i,j)来表 示，例如 A(3,2)就表示矩阵 A 的第三行第二列的元素. 可 以单独引用某个矩阵元素，并且对其直接赋值或建立表达 式进行修改. *（三）一些特殊矩阵的生成指令 建立 1 矩阵使用 ones 命令，这种矩阵的元素全部都是 1. 建立 0 矩阵使用 zeros 命令，这种矩阵的元素全部都是 0. 单位矩阵的对角线元素全部是 1，而其他元素全部是 0， 建立它使用 eye 命令. ones(n) 建立一个 n×n 的 1 矩阵. ones(size(A)) 建立一个和矩阵A同样大小的1矩阵. zeros(n) 建立一个 n×n 的 0 矩阵. zeros(size(A)) 建立一个和矩阵A同样大小的0矩阵. eye(n) 建立一个 n×n 的单位矩阵. eye(size(A)) 建立一个和矩阵A同样大小的单位矩 阵. 14 （四）数组与数组运算 对于一个矩阵 A，Matlab 既可以进行矩阵运算，也可 以进行数组运算. 所谓矩阵运算，就是线性代数中规定的 矩阵所特有的运算，例如矩阵的乘法，求矩阵的逆，特征 值等等；而数组运算，则是把参与运算的二维矩阵看作二 维的数组，对矩阵中的对应元素进行相同的运算，运算结 果到一个新的矩阵，新矩阵和原矩阵中的元素一一对应. 这种数组运算使得对大量的数据进行相同的运算变得十分 方便快捷，我们首先掌握数组运算. 数组运算的一般运算符与单个数值运算的一般运算符 相似，相应的有 + - .* ./ .^ .\ ( ) 含义与数值运算的算符相同. 例：比较“*”和“.*”的不同. >> A=[1,2,3;4,5,6;7,8,9]; >> B=[1,1,1;2,2,2;3,3,3]; >> D=A.*B D = 1 2 3 8 10 12 21 24 27 >> d=A*B d = 14 14 14 32 32 32 50 50 50 可见，两种算符得出的结果是不一样的，写表达式时 一定要分清要做哪种运算. 现在我们所做的通常是数组运 算. 参与运算的矩阵的维数必须相同，否则则会报错. 例 如： >> A=[1,2,3;4,5,6] A = 1 2 3 4 5 6 >> B=[1,2;3,4;5,6] B = 1 2 15 3 4 5 6 >> A.*B ??? Error using ==> .* Matrix dimensions must agree. 在Matlab中预定义的数学函数是基于矩阵的数组运算. 也就是说函数f对于矩阵A，有f (A)ij=f (Aij)，矩阵的维数没 有改变，例如： A = 1 2 3 4 5 6 >> sin(A).*cos(A).^2 ans = 0.2456 0.1575 0.1383 -0.3233 -0.0772 -0.2576 §5 编程 Matlab 输入命令的方式有两种，一种就是像前文所例 举的在命令窗口中直接输入简单的语句，这种方式适应于 命令比较简单、且处理的问题没有普遍应用性、差错处理 比较简单的场合. 但是在进行大量重复性的计算时，或者 语句结构比较复杂需要进行流程控制时，这种方式就显得 不太灵活，如果输入错误则很不好改正. 因此需要借助另 一种工作方式：M 文件的编程工作方式. M 文件是一个简 单的文本文件，语法比一般的高级语言都要简单，程序容 易调试，交互性强；而且可以像一般文本文件那样在任何 文本编辑器中进行编辑、存储、修改和读取（输入时用英 文,不可用中文标点,注意空格）. 前面所举的所有示例都可 以编辑为 M 文件. （一）创建并编辑 M 文件 尽管 M 文件可以用任何文本编辑器（例如记事本，写 字板等程序）创建和编辑，Matlab 还是提供了一个方便实 用的 M 文件编辑器，利用这个编辑器，用户可以完 16 成程序的创建、编辑、调试、存储和运行等工作. 在 Matlab 命令窗口中输入“edit”并回车，就调出了如 图 8 所示的 M 文件编辑器（编辑窗口）. 图 8 M 文件编辑器 图 9 修改参数λ后的运行结果 在 M 文 件 编 辑 器 的 菜 单 栏 里 选 择 “FileÆNewÆM-file”命令，或单击工具栏上的 按钮， 可以创建一个新的 M 文件；选择“FileÆOpen”命令，或 者单击工具栏上的 按钮，可以打开工作路径上的已有的 M 文件；当编辑完 M 文件后，在菜单栏里选择 “FileÆSave”命令或者单击工具栏上的 按钮，可以保 存文件. 此外，在 Matlab 主窗口中的菜单栏和工具栏上， 也有相应的命令和按钮来创建和编辑 M 文件. 例 1：将§2 中的 Logistic 叠代图的语句编辑成 M 文件并运行. 启动 M 文件编辑器，创建一个新的 M 文件，在编辑窗口中键 入下列语句： x0=abs(sin(randn)); lamda=0.8; xn=x0; for n=1:150 xn=lamda*xn*(1-xn); plot(n,xn, '.'); hold on end 输入完毕后保存文件为 Logistic.m，放在当前的工作目录下. 在 Matlab 的命令窗口中输入“Logistic”并回车，则运行 17 此程序，显示的结果与§2 节中相同. 用户随时可以更改 M 文件内容 并重新运行，例如将“lamda=0.8”改为“lamda=3.5”， 保存文件并再 次运行此程序，则所得结果如图 9 所示. 读者可以将前面所举的例子改为以 M 文件的形式来 完成. （二）脚本文件与函数文件 M 文件有两种形式，即脚本文件（也称命令文件）和 函数文件. 脚本文件就是前面所举例中的 M 文件形式，是 命令和语句的简单叠加，Matlab 会自动按顺序执行文件中 的命令. 函数文件主要用以解决参量传递和函数调用的问 题，它的第一句以 function 语句为引导. 需要注意的是命令文件在运行过程中可以调用Matlab 工作域内所有的数据，而且所产生的所有变量均为全局变 量. 也就是说，这些变量一旦生成，就一直保存在内存空 间中，直到用户执行“clear”命令或退出 Matlab 时才会被清 除. 而函数文件中所有变量除特殊声明外，均为局部变量， 即只在调用函数时存在，一旦函数调用结束，这些变量立 即被清除. 按下“F5”键或单击工具栏上的 按钮，可以运行该程 序. 例 2：将 Logistic 方程以函数文件形式完成. 在 M 文件编辑器中输入以下语句： function f=logistic(lamda,N) x0=abs(sin(randn)); xn=x0; for n=1:N xn=lamda*xn*(1-xn); plot(n,xn, '.'); hold on end 保存文件为 logistic.m，在命令窗口中输入 logistic(3.5,150)，所 得结果与图 9 相同. 18 在本例中，第一行语句必须以 function 开头，f 是函数的返回 值（可以是空值），logistic 是函数名，lamda 和 N 是参量，在调用此 函数文件时由用户给出. 之后的语句与前面例子相同，只是不再在 程序中对 lamda 和循环次数 N 赋值了. 调用函数文件的命令是 logsitic(3.5,150)，即向函数传递两个参 量（λ和 N），表示当 λ=3.5 时叠代 150 次，调用函数得出结果. 写 成函数文件的好处是不需要对 M 文件本身进行修改，只需在调用函 数时就可以随时改变参量的值，得到不同情况下的结果，十分方便. 例如在命令窗口中输入 logistic(3.2,200)，就是求 λ=3.2 时叠代 200 次的结果. 需要注意的是要保存函数文件的文件名与首行语句 中的函数名要一致，才能保证成功调用函数. （三）流程控制语句 对于比较复杂的任务，往往需要设计一定的算法，运 用流程控制的语句才能完成，因此有必要熟悉 Matlab 中的 流程控制语句和一些算法知识. 程序一般分为顺序结构、 循环结构、分支结构 3 种，在理论上，只要有以上三种结 构就可以构造功能强大的程序. 下面只介绍顺序结构 for 循环结构. （1）顺序结构 顺序结构就是依次顺序执行程序的各条语句. 语句在 程序文件中的物理位置就反映了程序的执行顺序. 例如绘 制等量异号点电荷电场的电势分布的语句就是典型的顺序 结构： [x,y]=meshgrid(-5:0.2:5,-4:0.2:4); %建立数据网格 z=1./sqrt((x-2).^2+y.^2+0.01)-1./sqrt((x+2).^2+y.^2+0.01); %电势的表达式 mesh(x,y,z) %三维曲面绘图 （2）for 循环结构 循环是计算机解决问题的主要手段，许许多多实际问 题大都包含有规律性的重复计算和对某些语句的重复执行. for 循环将循环体中的语句重复执行给定的次数，循 19 环的次数一般情况下是已知的，除非用其它的语句将循环 提前结束. 在前面介绍 Logistic 叠代的例子时就已经接触 到 for 循环了，其结构为： for i=表达式 可执行语句 1 …… 可执行语句 n end 表达式是一个数组，可以为 m:s:n，m,s,n 一般情况下 为整数，也可以取小数. 此数组中的元素被逐一赋值给 i， 对于每个 i 的不同取值都要执行一次循环体内的语句. 表 达式也可为 m:n，此时默认的步长为 1. for 循环的循环体中，可以多次嵌套 for 和其他的结构 体，这在一些情况下是非常有用的，大大扩展了 for 循环 的用途. 例 3：利用 for 循环求 1!+2!+3!+…+20!的值. （lt53.m） sum=0; %先将所求的值预定义为零 for i=1:20 %以 i 为循环变量，循环 20 次 prd=1; %每次循环中都设一个临时变量 prd，初始值为 1 for k=1:i %以 k 为循环变量，从 1 至 i 完成子循环 prd=prd*k; %计算 k 的阶乘 end %结束内层循环 sum=sum+prd; %把计算出的 k 的阶乘加入变量 sum 中 end %结束外层循环 sum %求值 sum = 2.5613e+018 循环的嵌套类似于时钟的运行，外层循环变量 i 是时 针，内层循环变量 k 是分针，当分针 k 完成一轮循环时， 时针 i 只完成了一步循环. *（3）其他指令 还有其它一些用于流程控制的命令，它们的名称和功 能如下表所示. pause 使程序运行暂停，按任意键恢复运行，而 pause(n)则是暂停 n 秒. 20 disp(‘…’) 在屏幕上显示引号中的内容 input 将用户从键盘输入的数值、字符串或表达式赋 予指定的变量. §6 普通物理学中常用的几种计算 本节介绍一些在普通物理中经常用到的计算指令和函 数. （一）解方程 （1）解一元函数方程 解一元函数方程 f(x)=0，即利用指令 fzero( )来求一元 函数 y=f (x)的零点. 语句格式及各项符号说明如下： x=fzero(fun,x0) x=fzero(fun,x,options) x=fzero(fun,x,options,p1,p2,…) [x,fval]=fzero(…) x 找到的零点，函数值在此会改变符号，没有零点则返 回 NaN. fval 零点对应的函数值（由于数值解法有一定的精确度， 所以往往零点对应的函数值不是恰好为零，而是一个 小量）. fun 用 M 文件或指令建立的单变量函数，其函数值必须是 实数. x0 猜测的初始值或搜寻零点的区间[x0(1),x0(2)]，必须是 实数. 如果是区间，要求函数值在 x0(1)与 x0(2)符号 相反. 如果是猜测值，则 fzero 将寻找一个包含 x0， 同时函数值符号发生改变的区间，如果寻找中遇到了 Inf，NaN 或复数，则寻找终止. 21 options 控制算法的优化选项，可选 display 和 TolX 两项，[] 是默认值. p1,p2 传递给函数的参数. 例 1：求方程 的根. 0cossin3cossin 22 =−+ xxxxx 首先需要在 M 文件编辑器中建立函数文件： function y=qx(x) y=x.*sin(x)+cos(x.^2)-3*sin(x).^2.*cos(x); 在当前工作路径下保存文件为 qx.m. 在 Matlab 命令窗口中输 入“[x,y]=fzero('qx',-3)”并回车，所得结果如下： >> [x,y]=fzero('qx',-3) x = -2.9396 y = 1.3878e-017 值得注意的是，用这种方法寻找函数的零点只能在初 始猜测值附近找到最接近的一个零点，如果函数有多个零 点，则有遗漏的可能. 所以想要找到函数的所有零点，最 好使用函数绘图语句 fplot( )，结合图像判断搜寻范围. 在 命令窗口中输入 fplot('qx',[-9 9])，可画出此函数在区间[-9 9] 上的图像，如图 10 所示. 可以看到在此区间上的-6、-3、3、 6 处均存在零点，可以分别以这几个点为猜测值代入指令， 求出较为精确的解. 图 10 y =xsinx+cosx2-3sin2xcosx在[-9 9]上的图像 *（2）解线性方程组 这里作的是矩阵运算，还是不必去学矩阵，你只要模仿例题的 做法，就能解线性方程组. 22 对形式为 Ax=b的线性方程组，这里 A是 m×m 的矩 阵，b是已知矢量，x是未知矢量. 用 Matlab 求解时使用 矩阵左除的方法，即使用程序式 x=A\b. 例 2：利用 Matlab 求下方程组的数值解. （lt62.m） ⎪⎩ ⎪⎨ ⎧ =−+ =+− =−+ 923 14354 28762 zyx zyx zyx >> A=[2 6 -7;4 -5 3; 3 1 -2]; %建立系数矩阵 >>b=[28;14;9]; %建立已知矢量 >>x=A\b %求解 x = -6.6471 -24.0000 -26.4706 本例程序中的 x 是未知数 x、y、z 所组成的解的矢量， 它的三个元素 x(1)、x(2)、x(3)分别代表未知数 x、y，z. （二）差分、梯度、积分 本小节介绍数值计算差分、梯度和积分的相关指令. （1）差分 给出矢量 x（数组），Matlab 提供了一个近似计算微分 的函数 diff，它计算数组中元素间的差分. diff 的用法为： diff(x) 其值为[x(2)-x(1);x(3)-x(2);…;x(n)-x(n-1)] 对函数 y=f(x)，计算导数 dy/dx 则用 diff(y)./diff(x)来进 行. 例 3：求 在[-π,π]上的导数.（lt63.m） xy sin= x=-pi:0.01*pi:pi; %在[-π,π]上每隔 0.01π取一个值 y=sin(x); %建立函数关系 dy=diff(y)./diff(x); %利用对 y 和 x 的差分计算各点上的导数 plot(x,y); %画出原函数曲线 hold on %图形控制，使绘制新图时不抹掉旧图 x(length(x)-1)=[]; %将 x 的最后一个元素去掉 plot(x,dy, 'r.') %绘制 x 和的 y 的函数关系，即导函数曲线 23 最终结果如图 11 所示，蓝线表示 y=sinx，红点为导函数曲线. 当 x 在区间上取值越密集时，用差分的办法求出的导数越精确 的. 如果自变量取值间隔太大，用有限差分近似计算微分会导致很 差的结果. 由于 diff 计算数组元素间的差分，输出的数组比原数组少了一 个元素. 所以，画导函数曲线时，必须舍弃 x 数组中的一个元素. 图 11 利用 diff 求导数 图 12 绘制等势面和电场强度 （2）梯度 Matlab 用 gradient( )函数来求矩阵梯度，其语句格式 如下： [fx,fy]=gradient(f) [fx,fy]=gradient(f,h) [fx,fy]=gradient(f,hx,hy) [fx,fy,fz]=gradient(f) [fx,fy,fz]=gradient(f,hx,hy,hz) 各项符号含义为： fx 代表 df/dx，是 x（列）方向的偏导数 fy 代表 df/dy，是 y（行）方向的偏导数 fz 代表 df/dz，是 z（层）方向的偏导数 f 求梯度的矩阵 h 标量，是各个方向的步长. 不指定 h 的时候则取默 认值 1 hx,hy,hz 分别表示在 x,y,z 方向上的步长 当 f 是矢量时，gradient 可以代替 diff 函数来求 f 的导 数. 例 3 也可以改为：（lt631.m） x=-pi:0.01*pi:pi; 24 y=sin(x); dy=gradient(y,0.01*pi); %用 gradient 函数求导，0.01π 为步长 plot(x,y); hold on plot(x,dy,'r.') 所得结果与例 3 没有明显区别. 例 4：对§2 中的“等量异号点电荷电场中电势分布”的示例，求 出沿 x,y方向电势的梯度，并画出等势面和各点上电场的大小和方向. （lt64.m） [x,y]=meshgrid(-2:0.1:2,-2:0.1:2); %以 0.1 为步长建立平面数据网格 z=1./sqrt((x-1).^2+y.^2+0.01)... %写出电势表达式 -1./sqrt((x+1).^2+y.^2+0.01); [px,py]=gradient(z); %求电势在 x,y 方向的梯度即电场强度 contour(x,y,z,[-12,-8,-5,-3,-1,... %画出等势线 -0.5,-0.1,0.1,0.5,1,3,5,8,12]) hold on %作图控制 quiver(x,y,px,py,'k') %画出各点上电场的大小和方向 计算结果如图 12 所示. 本例中为了减少计算量，增加精确度，与先前的示例相比，计 算范围由原先的-5> sum=quad('example',0,10) %调用 quad 函数求[0,10]上的积分 sum = 51.0000 （三）解常微分方程 Matlab 解常微分方程组的能力很强而且很方便，对于 我们在普通物理学中遇到的大多数动力学方程都可以用命 令 ode45 求解. Matlab 只能解一阶的常微分方程组，高阶的常微分方 程需要转化成一阶方程组才能求解. 对于二阶常微分方程 ,首先需要化成显式形式 ，然 后令 ， 0),,,( =txxxF &&& ),,( txxfx &&& = xy =)1( xy &=)2( ，则二阶常微分方程化为两个一 阶常微分方程组成的方程组，从而使问题得到解决. 26 )),2(),1(( d )2(d )2( d )1(d tyyf t y y t y = = 例 6：质点在万有引力作用下的运动. 以万有引力的固定不动 的施力质点m0所在位置为坐标原点O, 建立直角坐标系Oxy，质点的 运动微分方程为 r r Gmmrm r&&r 3 0−= ，分量方程为： 212222 0 )( yx x yx Gmx ++−=&& , 212222 0 )( yx y yx Gmy ++−=&& 这两个方程都是二阶常微分方程，定义解矢量为 y，令 xy =)1( , xy &=)2( , , yy =)3( yy &=)4( 可将方程组化为： )2( d )1(d y t y = , ( ) 2322 0 )3()1( )1( d )2(d yy yGm t y + ⋅−= )4( d )3(d y t y = , ( ) 2322 0 )3()1( )3( d )4(d yy yGm t y + ⋅−= 编写微分方程组函数文件 yxlcfun.m： function ydot=yxlcfun(t,y,flag,p)%函数首行，p为参量Gm0 ydot=[y(2); %建立微分方程组 p*y(1)/sqrt(y(1).^2+y(3).^2).^3; y(4); p*y(3)/sqrt(y(1).^2+y(3).^2).^3]; 解微分方程的主程序 yxlc.m： p=-1; %取Gm0=1 y0=[-10 0.2 6 0.2;-25 0.5 5 0;-25 0.8 6 0]; %三组不同的初始条件 plot(0,0,'*r') %画出 O 点 for i=1:3 %分别以不同初始条件解 3 次方程 [t,y]=ode45('yxlcfun',[0:0.1:300],y0(i,:),[ ],p); 27 hold on axis([-25 25 -20 20]); %指定坐标范围 comet(y(:,1),y(:,3)) %绘出质点运动轨迹(x,y) end %结束循环 解出的结果如图 13 所示. 由上面例子，我们初步了解了 Matlab 解常微分方程的 一般过程，首先是建立微分方程函数文件，文件的格式如 下： fuction ydot=filename(t,y, p1,p2) %t,y 是积分区间和解矩阵 %p1,p2 是参数 ydot=[关于 t,y 的表达式]; %ydot 表示 dy/dt 下面介绍 ode45 命令的用法，ode45 的一般调用格式 为： [T,Y]=ode45(‘fun’,tspan,y0,options,p1,p2,…) 其中含义为： fun 求解的微分方程函数名 tspan 单调递增（减）的积分区间[t0:tstep:tfinal] y0 初始条件矢量 options 用 odeset 建立的优化选项，一般都用默认值，为 空矢量“[ ]” p1,p2 传递给 fun 函数的参数 T,Y T 是输出的时间列矢量，矩阵 Y 的每一个列矢量 是解的一个分量 各个项在命令中的位置和顺序不能颠倒，否则程序就 会出错. 图 13 万有引力场中质点运动轨迹 图 14 图形窗口 28 §7 作图和动画 无论在科研还是数学上，计算的结果都希望能做到可 视化，而 Matlab 完善与强大的作图功能，正好满足了这种 需求. 本节介绍最基本的一些绘图指令和动画指令. （一）Matlab 的图形窗口 当使用绘图语句时，Matlab就自动打开一个图形窗口； 如果已经有图形窗口存在，作图命令便会使用已存在的图 形窗口. 如果使用命令 figure，就会打开一个新的图形窗口. 每个图形窗口的标题栏都会有一个编号 n，打开第 n 个图 形窗口的指令是 figure(n)；在已有图形上继续作图的指令 是 hold on；取消这种功能的指令是 hold off. 图形窗口如图 14 所示. 在图中，上端是菜单栏，选择 tools菜单中Show Toolbar 命令可以显示工具栏的图标. 各图标的含义分别是新建、 打开、保存和打印文件，编辑图形、加注文字、画箭头和 画直线，放大、缩小和旋转图形. 在按下工具栏上的“鼠标箭头”后，可以对当前的图形 对象的各种属性（颜色、线宽、坐标格式、标志等）加以 编辑. File 菜单下的“选项”（Preferences）可以改变图形窗 口的各种功能，如数据格式、字体，图形存储的格式等. （二）二维图形 （1）plot 语句 二维图形绘图命令中最基本的指令就是 plot，它有不 同的形式，与输入的量有关. 对于矢量 y，plot(y)产生一个 折线图，纵轴为 y 的元素，横轴为 y 的元素指标. 如果输 入两个矢量 x,y，则 plot(x,y)产生的是 y 相对于 x 的图形， 如前面列举的绘制振动曲线的例子. 29 一个 plot 命令可以同时画多个图形. Matlab 会自动用 不同颜色区分每组数据的图形，颜色是预先设置好的，但 可以改动. 指定颜色、线型和数据点标志的 plot 命令格式如下： plot(x1,y1,’color style marker’,x2,y2,’color style marker’,…) 例如 plot(x,y,’y:+’)画出的是黄色点状线，在每个数据 点有加号标志. 如果指定了数据点而不指定线型，则只在 各个数据点画出标志符号. 颜色、线型和数据点标志的可 选项有： 颜色 青’c’，洋红’m’，黄’y’，黑’k’，红’r’，绿’g’，蓝’b’，白’w’ 线型 实线’-’，短划线'--'，点线’:’，点划线’-.’ 标志 点’.’，加号’+’，圆圈’o’，星号’*’，叉号’×’，四方形’s’， 钻石形’d’，五角形’p’，六角形’h’，三角形（向上’v’， 向下’^’，向左’<’，向右’>’） 例 1：在同一窗口中画多个函数图像. （lt71.m） x=0.01*pi:0.01*pi:pi; %建立自变量数组 y1=(2*x-pi).^2-2; %各个函数表达式 y2=sqrt(abs(3*x)); y3=2*sin(2*x); y4=2*cos(2*x); y5=1./(10*x); figure; %开启图形窗口 plot(x,y1,x,y2,'g--.',x,y3,'m:*',x,y4,'k.',x,y5,'b-.') %绘制函数图像 运行结果如图 15 所示. （2）fplot 语句 当已知函数关系的时候，采用 fplot 语句可以更加快速 精确地绘制出指定区间上的函数图像. fplot 绘图的数据点 是自适应产生的，即在函数变化小的地方取较少的点，在 函数变化剧烈的地方取较多的点. 如 fplot(‘sin(1./x)’,[0.01 0.1],1e-3); 直接画出了在[0.01 0.1]区间上函数sin(1/x)的曲线，其 误差为 10-3，画出图形如图 16 所示. 30 图 15 用 plot 画出的函数图形 图 16 fplot 绘制的函数图像 读者可以尝试将例 1 的各函数用 fplot 语句绘制出来， fplot 的格