《数值运算(微分方程求解与函数优化).ppt》由会员分享，可在线阅读，更多相关《数值运算(微分方程求解与函数优化).ppt（81页珍藏版）》请在三一办公上搜索。
1、第二讲 MATLAB的数值计算,matlab 具有出色的数值计算能力，占据世界上数值计算软件的主导地位,数值运算的功能,创建矩阵矩阵运算多项式运算线性方程组数值统计线性插值函数优化微分方程的数值解,一、命令行的基本操作,创建矩阵的方法直接输入法规则：矩阵元素必须用 括住 矩阵元素必须用逗号或空格分隔 在 内矩阵的行与行之间必须 用分号分隔,a=1;b=2;c=3;x=5 b c;a*b a+c c/b x=5.000 2.000 3.000 2.000 4.000 1.500y=2,4,5;3 6 8 y=2 4 5 3 6 8,矩阵元素可以是任何matlab表达式，可以是实数，也可以是复数，
2、复数可用特殊函数I，j 输入。大的矩阵可以用分行输入，回车键代表分号。a=1 2 3;4 5 6 x=2 pi/2;sqrt(3)3+5i,矩阵元素,符号的作用,逗号和分号的作用 逗号和分号可作为指令间的分隔符，matlab允许多条语句在同一行出现。分号如果出现在指令后，屏幕上将不显示结果。,注意：只要是赋过值的变量，不管是否在屏幕上显示过，都存储在工作空间中，以后可随时显示或调用。变量名尽可能不要重复，否则会覆盖。当一个指令或矩阵太长时，可用续行,冒号的作用 用于生成等间隔的向量，默认间隔为1。用于选出矩阵指定行、列及元素。循环语句,2.用matlab函数创建矩阵,空阵 matlab允许输入
3、空阵，当一项操作无结果时，返回空阵。rand 随机矩阵eye 单位矩阵zeros 全部元素都为0的矩阵ones 全部元素都为1的矩阵diag 产生对角矩阵,例 eye(2，3)zeros(2，3)ans=ans=1 0 0 0 0 0 0 1 0 0 0 0 ones(2，3)ans=1 1 1 1 1 1 V=5 7 2；A=diag(V)A=5 0 0 0 7 0 0 0 2,例 eye(2)ans=1 0 0 1 zeros(2)ans=0 0 0 0 ones(2)ans=1 1 1 1,例 在区间20,50内均匀分布的5阶随机矩阵。命令如下：x=20+(50-20)*rand(5)此
4、外，常用的函数还有reshape(A,m,n)，它在矩阵总元素保持不变的前提下，将矩阵A重新排成mn的二维矩阵。,也可用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)、对角矩阵、范德蒙等矩阵的创建，就不一一介绍了。注意：matlab严格区分大小写字母，因此a与A是两个不同的变量。matlab函数名必须小写。,3.用m文件创建矩阵,对于比较大且比较复杂的矩阵，可以为它专门建立一
5、个M文件。下面通过一个简单例子来说明如何利用M文件创建矩阵。例 利用M文件建立MYMAT矩阵。,(1)启动有关编辑程序或Matlab文本编辑器，并输入待建矩阵。(2)把输入的内容以纯文本方式存盘(设文件名为mymatrix.m)。(3)在Matlab命令窗口中输入mymatrix，即运行该M文件，就会自动建立一个名为MYMAT的矩阵，可供以后使用。,4.用冒号表达式创建矩阵,利用冒号表达式可以线性等间距地建立一个向量来创建矩阵 一般格式是：e1:e2:e3其中e1为初始值，e2为步长，e3为终止值。或者为：（start:step:end）例 a=1:2:10 a=1 3 5 7 9,5.矩阵的
6、修改,直接修改 可用键找到所要修改的矩阵，用键移动到要修改的矩阵元素上即可修改。指令修改 可以用A(,)=来修改。,例如a=1 2 0;3 0 5;7 8 9a=1 2 0 3 0 5 7 8 9a(3,3)=0a=1 2 0 3 0 5 7 8 0,把Matlab工作空间中一些有用的数据长久保存下来的方法是生成mat数据文件。save 将工作空间中所有的变量存到matlab.mat文件中。,二、数据的保存与获取,默认文件名,save data将工作空间中所有的变量存到data.mat文件中。save data a b 将工作空间中a和b变量存到data.mat文件中。下次运行Matlab时即
7、可用load指令调用已生成的mat文件。,load load data load data a b mat文件是标准的二进制文件，还可以ASCII码形式保存。,即可恢复保存过的所有变量,矩阵加、减（,）运算规则：相加、减的两矩阵必须有相同的行和列两矩阵对应元素相加减。允许参与运算的两矩阵之一是标量。标量与矩阵的所有元素分别进行加减操作。,三、矩阵运算,2.矩阵乘（）运算规则：A矩阵的列数必须等于B矩阵的行数标量可与任何矩阵相乘a=1 2 3;4 5 6;7 8 0;b=1;2;3;c=a*bc=14 32 23,d=-1;0;2;f=pi*df=-3.1416 0 6.2832 矩阵除的运算在
8、线性代数中没有，有矩阵逆的运算，在matlab中有两种矩阵除运算。,两种除法：和/，分别表示左除和右除。如果A矩阵是非奇异方阵，则AB和B/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;a2 ans=30 36 42
9、66 81 96 102 126 150,当一个方阵有复数特征值或负实特征值时，非整数幂是复数阵。,a0.5 ans=0.4498+0.7623i 0.5526+0.2068i 0.6555-0.3487i 1.0185+0.0842i 1.2515+0.0228i 1.4844-0.0385i 1.5873-0.5940i 1.9503-0.1611i 2.3134+0.2717i,inv 矩阵求逆det 行列式的值eig 矩阵的特征值diag 对角矩阵 矩阵转置sqrt 矩阵开方,4.矩阵的其它运算,5.矩阵的范数,矩阵范数的函数为：(1)norm(V)或norm(V,2)：计算矩阵V的
10、2范数。(2)norm(V,1)：计算矩阵V的1范数。(3)norm(V,inf)：计算矩阵V的范数。,6.矩阵的一些特殊操作,矩阵的变维 a=1:12;b=reshape(a,3,4)c=zeros(3,4);c(:)=a(:)矩阵的变向 rot90:旋转;fliplr:左右翻;flipud:上下翻矩阵的抽取 diag:抽取主对角线;(对于非方阵的情况?)tril:抽取主下三角;triu:抽取主上三角，然后其余补零元素矩阵的扩展,关系运算,关系运算符的运算法则,当两个比较量是标量时，直接比较两数的大小。若关系成立，关系表达式结果为1，否则为0。当参与比较的量是两个维数相同的矩阵时，比较是对两
11、矩阵相同位置的元素按标量关系运算规则逐个进行，并给出元素比较结果。最终的关系运算的结果是一个维数与原矩阵相同的矩阵，它的元素由0或1组成。,(3)当参与比较的一个是标量，而另一个是矩阵时，则把标量与矩阵的每一个元素按标量关系运算规则逐个比较，并给出元素比较结果。最终的关系运算的结果是一个维数与原矩阵相同的矩阵，它的元素由0或1组成。注意：其书写方法与数学中的不等式符号不尽相同。,数组运算指元素对元素的算术运算，与通常意义上的由符号表示的线性代数矩阵运算不同。数组加减(.+,.-)a.+b a.-b,7.矩阵的数组运算,对应元素相加减（与矩阵加减等效）,2.数组乘除(，./，.）ab a，b两数
12、组必须有相同的行 和列两数组相应元素相乘。a=1 2 3;4 5 6;7 8 9;b=2 4 6;1 3 5;7 9 10;a.*bans=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 9 10;a*bans=25 37 46 55 85 109 85 133 172,a./b=b.aa.b=b./aa./b=b.a 都是a的元素被b的对应元 素除，“/”是斜杠a.b=b./a 都是b的元素被a的对应元 素除，“”是反斜杠例:a=1 2 3;b=4 5 6;c1=a.b;c2=b./ac1=4.0000 2.5000 2
13、.0000c2=4.0000 2.5000 2.0000,给出a,b对应元素间的商.,3.数组乘方(.)元素对元素的幂例:a=1 2 3;b=4 5 6;z=a.2z=1.00 4.00 9.00z=a.bz=1.00 32.00 729.00(1.4 2.5 3.6),matlab语言把多项式表达成一个行向量，该向量中的元素是按多项式降幂排列的。f(x)=a0 xn+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;
14、p=poly(a)p=1.00-6.00-72.00-27.00p是多项式p(x)=x3-6x2-72x-27的系数matlab描述方法，我们可用：p1=poly2str(p,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.00r=roots(p)-求由p构成的多项式的根r=12.12-5.73 显然 r是矩阵a的特征值-0.39,当然我们可用poly令其返回多项式形式(这是poly的第二个功能)p2=pol
15、y(r)p2=1.00-6.00-72.00-27.00matlab规定多项式系数向量用行向量表示，一组根用列向量表示。,P=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多项式乘运算(向量卷积),例
16、:a(x)=x2+2x+3;b(x)=4x2+5x+6;c=(x2+2x+3)(4x2+5x+6)a=1 2 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.00p=poly2str(c,x)其中x表示自变量 p=4 x4+13 x3+28 x2+27 x+18,4.deconv多项式除运算(解卷积),a=1 2 3;c=4.00 13.00 28.00 27.00 18.00d=deconv(c,a)d=4.00 5.00 6.00,它们之间的关系为:c=conv(a,d)+r,5.多项式导数或微分,m
17、atlab提供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)ans=x4+2 x3+3 x2+4 x+5b=polyder(a)b=4 6 6 4poly2str(b,x)ans=4 x3+6 x2+6 x+4,6.多项式的积分,matlab提供polyint函数计算多项式的积分。命令格式：polyint(p,k):求多项式p的积分，设积分的常数项为k，polyint(p)默
18、认k=0 例：a=1 2 3 4 5;poly2str(a,x)ans=x4+2 x3+3 x2+4 x+5b=polyint(a,8)b=0.2 0.5 1.0 2.0 5.0 8.0poly2str(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
19、矩阵逆两种解:x=inv(a)b 采用求逆运算解方程 x=ab 采用左除运算解方程 注：若a为奇异的，则Matlab适当给出警告信息或者给出结果为inf。,方程ax=ba=1 2;2 3;b=8;13;x=inv(a)*b x=ab x=x=2.00 2.00 3.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=2 3x1
20、+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 0 0.00,=,a x=b,3.欠定方程组的解,当方程数少于未知量个数时,即不定情况,有无穷多个解存在。matlab可求出两个解：用除法求的解x是具有最多零元素的解是具有最小长度或范数的解，这个解是基于伪逆pinv求得的。,x1+2x2+3x3=1 2x1+3x2+4x3=2a=1 2 3;2 3 4;b=1;2;x=ab x=pinv(a)b x=x=1.00 0.83 0 0.33 0-0.17,a x=b,六、微分方程求解,微分方程求解的仿真算法有多种，常
21、用的有Euler（欧拉法）、Runge Kutta(龙 格-库塔法。Euler法称一步法，用于一阶微分方程,当给定仿真步长时：所以 yn+1=yn+hf(xn,yn)n=0,1,2 y(x0)=y0,Runge Kutta法龙格-库塔法：实际上取两点斜率的平均 斜率来计算的，其精度高于欧拉算法。龙格-库塔法：ode23 ode45,k1=hf(xn,yn)k2=hf(xn+h,yn+k1),例：x+(x2-1)x+x=0为方便令x1=x，x2=x分别对x1,x2求一阶导数，整理后写成一阶微分方程组形式 x1=x2 x2=x2(1-x12)-x1建立m文件解微分方程,建立m文件function
22、xdot=wf(t,x)xdot=zeros(2,1)xdot(1)=x(2)xdot(2)=x(2)*(1-x(1)2)-x(1)给定区间、初始值;求解微分方程t0=0;tf=20;x0=0 0.25;t,x=ode23(wf,t0,tf,x0)plot(t,x),figure(2),plot(x(:,1),x(:,2),命令格式:T,Y=ODE23(ODEFUN,TSPAN,Y0)建立m文件function dxdt=wf(t,x)dxdt=x(2);x(2)*(1-x(1)2)-x(1);求解微分方程t,x=ode23(wf,0 30,0 0.25);plot(t,x);figure(2
23、)plot(x(:,1),x(:,2),七、函数优化寻优函数：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处有最小值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对矩阵操作的规定：如果是向量
24、，则对数据整体操作；如果是矩阵，则对矩阵的列操作。max 各列最大值 mean 各列平均值sum 各列求和std 各列标准差var 各列方差sort 各列递增排序cumsum 元素累计和cumprod 元素累计积,八、数据分析,数据分析相关的函数位于目录：toolboxsmatlabdatafun下。Matlab对矩阵操作的规定：如果是向量，则对数据整体操作；如果是矩阵，则对矩阵的列操作。max 各列最大值 mean 各列平均值sum 各列求和std 各列标准差var 各列方差sort 各列递增排序cumsum 元素累计和cumprod 元素累计积,例：x=1 3 2 4,y=1 2 3 8;
25、5 6 7 4sort(x),sort(y),max(y)1 2 3 4 1 2 3 4 5 6 7 8 5 6 7 8,九、拟合与插值,多项式拟合 采用最小二乘法对给定的数据进行多项式拟合，最后给出多项式的系数。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.6915-87.1174 40.0070-0.9043xx=0:0.01:1;yy=poly
26、val(p,xx);plot(xx,yy,-b,x0,y0,or),曲线拟合图形用户接口,Matlab7.0提供了支持曲线拟合的图形用户接口。在“Figure”窗口“ToolsBasic Fitting”菜单中。为了使用该工具，先用待拟合的数据画图。x=0:0.2:10;y=0.25*x+20*sin(x);plot(x,y,ro);在复选框“Plot fits”中选择“cubic”。,2.插值插值的定义是对某些集合给定的数据点之间函数的估值方法。当不能很快地求出所需中间点的函数时，插值是一个非常有价值的工具，它可以在已知数据之间寻找估计值，常用到信号处理和图像处理中。Matlab提供了一维、
27、二维、三次样条等许多插值选择。,interp1一维插值 interp2二维插值interp3三维插值spline三次样条插值griddata 栅格数据插值利用已知点确定未知点粗糙 精确集合大的 简化的,一维插值就是对一维函数y=f(x)进行插值。yi=interp1(x,y,xi,method),x必须是向量，y可以是向量也可以是矩阵。(x,y)代表的是已知数据。这时，xi代表需要估计值的位置，yi表示插值后的估计值。method用于指定插值的方法：1.method=nearest，在已知数据的最临近点设置插值点，对插值点的数进行四舍五入。对超出范围的点返回一个NaN。此方法是最快的插值方法，
28、但数据平滑方面最差，其得到的数据是不连续的。,2.method=linear，采用直线连接相邻的两点，即线性插值，是此函数的缺省默认方法。执行速度比最临近插值稍慢，数据平滑要由于临近插值，且数据是连续的。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,nea
29、rest)yi3=interp1(x,y,xi,spline)yi4=interp1(x,y,xi,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同样可以采用最临近插值、双线性插值、三次样条插值。,小 结 本节介绍了matlab语言的数值运算功能，通过学习应该掌握：如何创建矩阵、修改矩阵符号的用法矩阵及数组运算多项式运算线性方程组与微分运算,
