1、第二讲 MATLAB的数值计算, matlab 具有出色的数值计算能力，占据世界上数值计算软件的主导地位,数值运算的功能,创建矩阵 矩阵运算 多项式运算 线性方程组 数值统计 线性插值 函数优化 微分方程的数值解,一、命令行的基本操作,创建矩阵的方法 直接输入法,a=1; b=2; c=3; x=5 b c; a*b a+c c/bx=5.000 2.000 3.0002.000 4.000 1.500 y=2, 4, 5;3 6 8 y= 2 4 53 6 8,注意：只要是赋过值的变量，不管是否在屏幕上显示过，都存储在工作空间中，以后可随时显示或调用。变量名尽可能不要重复，否则会覆盖 。当一
2、个指令或矩阵太长时，可用续行,2.用matlab函数创建矩阵,空阵 matlab允许输入空阵，当一项操作无结果时，返回空阵。 rand 随机矩阵 eye 单位矩阵 zeros 全部元素都为0的矩阵 ones 全部元素都为1的矩阵 diag 产生对角矩阵,也可用linspace函数产生行向量。其调用格式为： linspace(a, b, n) 其中a和b是生成向量的第一个和最后一个元素，n是元素总数。 例 a=linspace(1 , 10 , 10)a=1 2 3 4 5 6 7 8 9 10,还有伴随矩阵、稀疏矩阵、魔方矩阵(magic)、对角矩阵、范德蒙等矩阵的创建，就不一一介绍了。 注意
3、：matlab严格区分大小写字母，因此a与A是两个不同的变量。matlab函数名必须小写。,把Matlab工作空间中一些有用的数据长久保存下来的方法是生成mat数据文件。save 将工作空间中所有的变量存到matlab.mat文件中。,二、数据的保存与获取,默认文件名,save data将工作空间中所有的变量存到data.mat文件中。 save data a b 将工作空间中a和b变量存到data.mat文件中。下次运行Matlab时即可用load指令调用已生成的mat文件。,load load data load data a b mat文件是标准的二进制文件，还可以ASCII码形式保存。
4、,即可恢复保存过的所有变量,矩阵加、减（,）运算 规则： 相加、减的两矩阵必须有相同的行和列两矩阵对应元素相加减。 允许参与运算的两矩阵之一是标量。标量与矩阵的所有元素分别进行加减操作。,三、矩阵运算,2. 矩阵乘（）运算 规则： A矩阵的列数必须等于B矩阵的行数 标量可与任何矩阵相乘 a=1 2 3;4 5 6;7 8 0;b=1;2;3;c=a*b c =143223,d=-1;0;2; f=pi*d f = -3.141606.2832矩阵除的运算在线性代数中没有，有矩阵逆的运算，在matlab中有两种矩阵除运算。,两种除法：和/，分别表示左除和右除。如果A矩阵是非奇异方阵，则AB和B/
5、A运算可以实现。 AB等效于A的逆左乘B矩阵，而B/A等效于A矩阵的逆右乘B矩阵。 对于矩阵来说，左除和右除表示两种不同的除数矩阵和被除数矩阵的关系。对于矩阵运算，一般ABB/A。,a p a 自乘p次幂,方阵,1的整数,3. 矩阵乘方 an,ap,pa,对于p的其它值,计算将涉及特征值 和特征向量，如果p是矩阵，a是标量 ap使用特征值和特征向量自乘到p次 幂；如a,p都是矩阵，ap则无意义。,a=1,2,3;4,5,6;7,8,9;a2ans =30 36 4266 81 96102 126 150,当一个方阵有复数特征值或负实特征值时，非整数幂是复数阵。,inv 矩阵求逆 det 行列式
6、的值 eig 矩阵的特征值 diag 对角矩阵 矩阵转置 sqrt 矩阵开方,4. 矩阵的其它运算,数组运算指元素对元素的算术运算， 与通常意义上的由符号表示的线性代数 矩阵运算不同。数组加减(.+,.-)a.+ba.- b,7. 矩阵的数组运算,对应元素相加减（与矩阵加减等效）,2. 数组乘除(，./，.） ab a，b两数组必须有相同的行 和列两数组相应元素相乘。 a=1 2 3;4 5 6;7 8 9; b=2 4 6;1 3 5;7 9 10; a.*b ans =2 8 18 4 15 30 49 72 90,a=1 2 3;4 5 6;7 8 9; b=2 4 6;1 3 5;7
7、9 10;a*b ans =25 37 46 55 85 109 85 133 172,a./b=b.a a.b=b./a a./b=b.a 都是a的元素被b的对应元 素除， “/”是斜杠 a.b=b./a 都是b的元素被a的对应元素除， “”是反斜杠 例: a=1 2 3;b=4 5 6; c1=a.b; c2=b./a c1 = 4.0000 2.5000 2.0000 c2 = 4.0000 2.5000 2.0000, 给出a,b对应元素间的商.,3. 数组乘方(.) 元素对元素的幂 例: a=1 2 3;b=4 5 6; z=a.2 z =1.00 4.00 9.00 z=a.b z
8、 =1.00 32.00 729.00(1 .4 2 .5 3 .6),matlab语言把多项式表达成一个行向量， 该向量中的元素是按多项式降幂排列的。f(x)=a0xn+a1xn-1+an-1x+an可用行向量 p=a0 a1 an-1 an 表示 poly 产生特征多项式系数向量 特征多项式一定是n+1维的 特征多项式第一个元素一定是1,四、 多项式运算,例:a=1 2 3;4 5 6;7 8 0; p=poly(a) p =1.00 -6.00 -72.00 -27.00 p是多项式p(x)=x3-6x2-72x-27的系数 matlab描述方法，我们可用： p1=poly2str(p,
9、x) 函数文件，显示 数学多项式的形式 p1 =x3 - 6 x2 - 72 x 27 注意：多项式中缺少的幂次用0补齐。,2.roots 求多项式的根,a=1 2 3;4 5 6;7 8 0;p=poly(a) p =1.00 -6.00 -72.00 -27.00 r=roots(p)-求由p构成的多项式的根 r = 12.12-5.73 显然 r是矩阵a的特征值-0.39,当然我们可用poly令其返回多项式形式(这是poly的第二个功能) p2=poly(r) p2 =1.00 -6.00 -72.00 -27.00 matlab规定多项式系数向量用行向量表示，一组根用列向量表示。,P=
10、poly(r),输入r是多项式所有根，返回值为代表多项式的行向量形式。 P=poly(A)，输入是N*N的方阵，返回值p是长度为N+1的行向量多项式，它是矩阵A的特征多项式，也就是说多项式p的根是矩阵A的特征值。,求根的另一种方法,str1=x3-6x2-72x-27; p1=str2poly(str1); r=roots(p1);注：str2poly 实现把一个字符串表示的多项式转换为一个行向量表示的多项式。poly2str 同理。,3.conv多项式乘运算(向量卷积),例:a(x)=x2+2x+3; b(x)=4x2+5x+6; c = (x2+2x+3)(4x2+5x+6) a=1 2
11、3;b=4 5 6; c=conv(a,b)或c=conv(1 2 3,4 5 6) c = 4.00 13.00 28.00 27.00 18.00 p=poly2str(c,x) 其中x表示自变量 p = 4 x4 + 13 x3 + 28 x2 + 27 x + 18,5.多项式导数或微分,matlab提供polyder函数计算多项式的导数。 命令格式： polyder(p): 求p的导数 polyder(a,b): 求多项式a,b乘积的导数 p,q=polyder(a,b): 求多项式a除以b的商的导数，并以p/q的格式表示。,例：a=1 2 3 4 5; poly2str(a,x)
12、ans = x4 + 2 x3 + 3 x2 + 4 x + 5 b=polyder(a) b = 4 6 6 4 poly2str(b,x) ans =4 x3 + 6 x2 + 6 x + 4,6.多项式的积分,matlab提供polyint函数计算多项式的积分。 命令格式： polyint(p,k): 求多项式p的积分，设积分的常数项为k， polyint(p) 默认k=0 例：a=1 2 3 4 5; poly2str(a,x) ans = x4 + 2 x3 + 3 x2 + 4 x + 5 b=polyint(a,8) b = 0.2 0.5 1.0 2.0 5.0 8.0 pol
13、y2str(b,x) ans =0.2 x5 + 0.5 x4 + x3 + 2 x2 + 5 x + 8,五、代数方程组求解,matlab中有两种除运算左除和右除。 对于方程ax=b，a 为anm矩阵，有三种情况： 当n=m时，此方程成为“恰定”方程 当nm时，此方程成为“超定”方程 当nm时，此方程成为“欠定”方程matlab定义的除运算可以很方便地解上 述三种方程,1.恰定方程组的解,方程ax=b (a为非奇异)x=a-1 b 矩阵逆 两种解: x=inv(a)b 采用求逆运算解方程 x=ab 采用左除运算解方程 注：若a为奇异的，则Matlab适当给出警告信息或者给出结果为inf。,方
14、程ax=b a=1 2;2 3;b=8;13; x=inv(a)*b x=abx = x = 2.00 2.003.00 3.00,=,a x = b,例: x1+2x2=8 2x1+3x2=13,2.超定方程组的解,方程的个数大于未知量个数时，方程一般无解。 方程解 (a a)x=a b x=(a a)-1 a b 求逆法 (也用到了最小二乘解的原理) x=ab matlab用最小二乘法找一个准确地基本解。,例: x1+2x2=1 2x1+3x2=23x1+4x2=3a=1 2;2 3;3 4;b=1;2;3;解1 x=ab 解2 x=inv(aa) a b x = x =1.00 1.00
15、0 0.00,=,a x = b,3.欠定方程组的解,当方程数少于未知量个数时,即不定 情况,有无穷多个解存在。 matlab可求出两个解： 用除法求的解x是具有最多零元素的解 是具有最小长度或范数的解，这个解是基于伪逆pinv求得的。,x1+2x2+3x3=12x1+3x2+4x3=2a=1 2 3;2 3 4;b=1;2;x=ab x=pinv(a)b x = x =1.00 0.830 0.330 -0.17,a x = b,六、微分方程求解,微分方程求解的仿真算法有多种，常用的有Euler（欧拉法）、Runge Kutta(龙格-库塔法。 Euler法称一步法，用于一阶微分方程,当给定
16、仿真步长时：所以 yn+1 = yn + hf (xn,yn) n=0,1,2y(x0)=y0,Runge Kutta法龙格-库塔法：实际上取两点斜率 的平均 斜率来计算的，其精度高 于欧拉算法 。 龙格-库塔法：ode23 ode45,k1=hf(xn,yn) k2=hf(xn+h,yn+k1),七、函数优化 寻优函数： fmin 单变量函数 fmins 多变量函数 constr 有约束条件,无约束条件,例1：f(x)=x2+3x+2在-5 5区间的最小值 f=fmin(x2+3*x+2,-5,5) 例2：f(x)=100(x2-x12)2+(a-x1)2在x1=a, x2=a2处有最小值
17、function f=xun(x,a) f=100*(x(2)-x(1).2).2+(a-x(1).2; x=fmins(xun,0,0, , ,sqrt(2),八、数据分析,数据分析相关的函数位于目录：toolboxsmatlabdatafun下。 Matlab对矩阵操作的规定：如果是向量，则对数据整体操作；如果是矩阵，则对矩阵的列操作。 max 各列最大值 mean 各列平均值 sum 各列求和 std 各列标准差 var 各列方差 sort 各列递增排序 cumsum 元素累计和 cumprod 元素累计积,例：x=1 3 2 4,y=1 2 3 8;5 6 7 4 sort(x),so
18、rt(y),max(y)1 2 3 41 2 3 45 6 7 85 6 7 8,九、计算方法拟合与插值,1、拟合,目的是在众多的样本点中进行拟合，找出满足样本点的分布。,在科学研究和工程技术问题中，常需要研究变量与变量之间的关系。它们的关系通常表现为：,1、确定性关系 如 I=VR (欧姆定律),2、非确定性关系,由于因素的复杂性或其它原因，变量之间找不出完全确定的关系，我们只好舍而求其次，寻求变量之间的近似关系,拟 合,2.拟合的基本原理,1. 拟合问题引例,拟 合 问 题 引 例 2,求血药浓度随时间的变化规律c(t).,作半对数坐标系下的图形,曲 线 拟 合 问 题 的 提 法,已知一
19、组（二维）数据，即平面上 n个点（xi,yi) i=1,n, 寻求一个函数（曲线）y=f(x), 使 f(x) 在某种准则下与所有数据点最为接近，即曲线拟合得最好。,y=f(x),i 为点（xi,yi) 与曲线 y=f(x) 的距离,曲线拟合问题最常用的解法线性最小二乘法的基本思路,第一步:先选定一组函数 r1(x), r2(x), rm(x), mn, 令f(x)=a1r1(x)+a2r2(x)+ +amrm(x) （1） 其中 a1,a2, am 为待定系数。,第二步: 确定a1,a2, am 的准则（最小二乘准则）： 使n个点（xi,yi) 与曲线 y=f(x) 的距离i 的平方和最小
20、。,记,问题归结为，求 a1,a2, am 使 J(a1,a2, am) 最小。,从给出的离散数据（如实验数据、统计数据等）中找出变量之间的近似关系问题称为数据拟合问题,解决这类问题的方法称为数据拟合法。 最常用的是最小二乘拟合法,采用最小二乘法对给定的数据进行多项式拟合，最后给出多项式的系数。p=polyfit(x,y,n)，采用n次多项式p来拟合数据x和y，从而使得y与p(x)最小均方差最小。,x0=0:0.1:1; y0=-.447 1.978 3.11 5.25 5.02 4.66 4.01 4.58 3.45 5.35 9.22; p=polyfit(x0,y0,3) p = 56.
21、6915 -87.1174 40.0070 -0.9043 xx=0:0.01:1;yy=polyval(p,xx); plot(xx,yy,-b,x0,y0,or),2、插值,在解决实际问题时，用来描述客观现象的函数通常是很复杂的，往往很难写出具体表达式，如通过实验或测量得到了函数 y=f(x) 在一系列互异点 处的值 ，需要构造一个简单函数 作为函数 y=f(x) 的近似表达式,这类问题称为插值问题。,一维插值的定义,返回,拉格朗日插值,分段线性插值,三次样条插值,二、插值的方法,Matlab提供了一维、二维、 三次样条等许多插值选择。,interp1一维插值 interp2二维插值 in
22、terp3三维插值 spline三次样条插值 griddata 栅格数据插值,一维插值就是对一维函数y=f(x)进行插值。 yi=interp1(x,y,xi,method),x必须是向量，y可以是向量也可以是矩阵。(x,y)代表的是已知数据。这时，xi代表需要估计值的位置，yi表示插值后的估计值。method用于指定插值的方法： 1.method=nearest，在已知数据的最临近点设置插值点，对插值点的数进行四舍五入。对超出范围的点返回一个NaN。此方法是最快的插值方法，但数据平滑方面最差，其得到的数据是不连续的。,2.method=linear，采用直线连接相邻的两点，即线性插值，是此函
23、数的缺省默认方法。执行速度比最临近插值稍慢，数据平滑要由于临近插值，且数据是连续的。 3. method=spline，采用三次样条函数来获得插值点。处理速度最慢，可以产生最光滑的结果。Matlab提供了一个样条插值工具箱，位于 toolboxsplines下。 4. method=pchip，采用分段三次埃尔米特多项式插值。,例：x=0:2*pi;y=sin(x);xi=0:0.1:8; yi1=interp1(x,y,xi, linear) yi2=interp1(x,y,xi, nearest) yi3=interp1(x,y,xi, spline) yi4=interp1(x,y,xi
24、, cubic) p=polyfit(x,y,3);yy=polyval(p,xi); subplot(3,2,1);plot(x,y,o); subplot(3,2,2);plot(x,y,o,xi,yy); subplot(3,2,3);plot(x,y,o,xi,yi1); subplot(3,2,4);plot(x,y,o,xi,yi2); subplot(3,2,5);plot(x,y,o,xi,yi3); subplot(3,2,6);plot(x,y,o,xi,yi4);,二维插值,二维插值主要应用于图像处理和数据的可视化，对双变量的函数z=f(x,y)进行插值。 zi=interp2(x,y,z,xi,yi,method),原始数据x,y,z决定插值函数z=f(x,y),返回值zi是(xi,yi)在函数f(x,y)上的值。 method同样可以采用最临近插值、双线性插值、三次样条插值。,
