数值运算(微分方程求解与函数优化).ppt

第二讲 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.500 y=2, 4, 5; 3 6 8 y= 2 4 5 3 6 8,矩阵元素可以是任何matlab表达式 ，可以是实数 ，也可以是复数，复数可用特殊函数I，j 输入。大的矩阵可以用分行输入，回车键代表分号。 a=1 2 3;4 5 6 x=2 pi/2; sqrt(3) 3+5i,矩阵元素,符号的作用,逗号和分号的作用 逗号和分号可作为指令间的分隔符，matlab允许多条语句在同一行出现。 分号如果出现在指令后，屏幕上将不显示结果。,注意：只要是赋过值的变量，不管是否在屏幕上显示过，都存储在工作空间中，以后可随时显示或调用。变量名尽可能不要重复，否则会覆盖 。 当一个指令或矩阵太长时，可用续行,冒号的作用 用于生成等间隔的向量，默认间隔为1。 用于选出矩阵指定行、列及元素。 循环语句,2.用matlab函数创建矩阵,空阵 matlab允许输入空阵，当一项操作无结果时，返回空阵。 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) 此外，常用的函数还有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文件创建矩阵,对于比较大且比较复杂的矩阵，可以为它专门建立一个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. 矩阵的修改, 直接修改 可用键找到所要修改的矩阵，用键移动到要修改的矩阵元素上即可修改。 指令修改 可以用A(,)= 来修改。,例如 a=1 2 0;3 0 5;7 8 9 a =1 2 0 3 0 5 7 8 9 a(3,3)=0 a =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时即可用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*b c =14 32 23,d=-1;0;2; f=pi*d f = -3.1416 0 6.2832 矩阵除的运算在线性代数中没有，有矩阵逆的运算，在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 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的 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。 当参与比较的量是两个维数相同的矩阵时，比较是对两矩阵相同位置的元素按标量关系运算规则逐个进行，并给出元素比较结果。最终的关系运算的结果是一个维数与原矩阵相同的矩阵，它的元素由0或1组成。,(3) 当参与比较的一个是标量，而另一个是矩阵时，则把标量与矩阵的每一个元素按标量关系运算规则逐个比较，并给出元素比较结果。最终的关系运算的结果是一个维数与原矩阵相同的矩阵，它的元素由0或1组成。 注意：其书写方法与数学中的不等式符号不尽相同。,数组运算指元素对元素的算术运算， 与通常意义上的由符号表示的线性代数 矩阵运算不同。 数组加减(.+,.-) a.+b a.- 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 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 = 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; p=poly(a) p =1.00 -6.00 -72.00 -27.00 p是多项式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.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=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 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,4.deconv多项式除运算(解卷积),a=1 2 3; c = 4.00 13.00 28.00 27.00 18.00 d=deconv(c,a) d =4.00 5.00 6.00,它们之间的关系为: c = conv(a,d)+r,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) 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 poly2str(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。,方程ax=b a=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+4x2=3 a=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=2 a=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,六、微分方程求解,微分方程求解的仿真算法有多种，常用的有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 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) 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对矩阵操作的规定：如果是向量，则对数据整体操作；如果是矩阵，则对矩阵的列操作。 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;5 6 7 4 sort(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.9043 xx=0:0.01:1;yy=polyval(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提供了一维、二维、 三次样条等许多插值选择。,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。此方法是最快的插值方法，但数据平滑方面最差，其得到的数据是不连续的。,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, nearest) 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语言的数值运算 功能，通过学习应该掌握： 如何创建矩阵、修改矩阵 符号的用法 矩阵及数组运算 多项式运算 线性方程组与微分运算,