数值运算微分方程求解与函数优化

关于数值运算微分方程求解与函数优化第1页，共81页，2022年，5月20日，23点23分，星期四一、命令行的基本操作创建矩阵的方法直接输入法规则：矩阵元素必须用[]括住矩阵元素必须用逗号或空格分隔在[]内矩阵的行与行之间必须用分号分隔

第2页，共81页，2022年，5月20日，23点23分，星期四》a=1;b=2;c=3;》x=[5bc;a*ba+cc/b]x=5.0002.0003.0002.0004.0001.500》y=[2,4,5;368]y=245368第3页，共81页，2022年，5月20日，23点23分，星期四

矩阵元素可以是任何matlab表达式，可以是实数，也可以是复数，复数可用特殊函数I，j输入。大的矩阵可以用分行输入，回车键代表分号。

a=[123;456]x=[2pi/2;sqrt(3)3+5i]

矩阵元素第4页，共81页，2022年，5月20日，23点23分，星期四符号的作用逗号和分号的作用

逗号和分号可作为指令间的分隔符，matlab允许多条语句在同一行出现。

分号如果出现在指令后，屏幕上将不显示结果。第5页，共81页，2022年，5月20日，23点23分，星期四注意：只要是赋过值的变量，不管是否在屏幕上显示过，都存储在工作空间中，以后可随时显示或调用。变量名尽可能不要重复，否则会覆盖。当一个指令或矩阵太长时，可用•••续行第6页，共81页，2022年，5月20日，23点23分，星期四冒号的作用

用于生成等间隔的向量，默认间隔为1。用于选出矩阵指定行、列及元素。循环语句第7页，共81页，2022年，5月20日，23点23分，星期四2.用matlab函数创建矩阵空阵[]—matlab允许输入空阵，当一项操作无结果时，返回空阵。rand——随机矩阵eye——单位矩阵zeros——全部元素都为0的矩阵ones——全部元素都为1的矩阵diag——产生对角矩阵第8页，共81页，2022年，5月20日，23点23分，星期四例

》eye(2，3)》zeros(2，3)ans=ans=

1

0

00

0

00

1

00

0

0》ones(2，3)ans=1

1

11

1

1》V=[572]；A=diag(V)A=5

0

00

7

00

0

2第9页，共81页，2022年，5月20日，23点23分，星期四例

》eye(2)

ans=

1

0

0

1》zeros(2)

ans=

0

0

0

0》ones(2)

ans=

1

1

1

1第10页，共81页，2022年，5月20日，23点23分，星期四例在区间[20,50]内均匀分布的5阶随机矩阵。命令如下：

x=20+(50-20)*rand(5)此外，常用的函数还有reshape(A,m,n)，它在矩阵总元素保持不变的前提下，将矩阵A重新排成m×n的二维矩阵。第11页，共81页，2022年，5月20日，23点23分，星期四也可用linspace函数产生行向量。其调用格式为：linspace(a,b,n)其中a和b是生成向量的第一个和最后一个元素，n是元素总数。例

》a=linspace(1,10,10)a=

12345678910第12页，共81页，2022年，5月20日，23点23分，星期四

还有伴随矩阵、稀疏矩阵、魔方矩阵(magic)、对角矩阵、范德蒙等矩阵的创建，就不一一介绍了。注意：matlab严格区分大小写字母，因此a与A是两个不同的变量。

matlab函数名必须小写。第13页，共81页，2022年，5月20日，23点23分，星期四3.用m文件创建矩阵对于比较大且比较复杂的矩阵，可以为它专门建立一个M文件。下面通过一个简单例子来说明如何利用M文件创建矩阵。例利用M文件建立MYMAT矩阵。第14页，共81页，2022年，5月20日，23点23分，星期四(1)启动有关编辑程序或Matlab文本编辑器，并输入待建矩阵。(2)把输入的内容以纯文本方式存盘(设文件名为mymatrix.m)。(3)在Matlab命令窗口中输入mymatrix，即运行该M文件，就会自动建立一个名为MYMAT的矩阵，可供以后使用。第15页，共81页，2022年，5月20日，23点23分，星期四4.用冒号表达式创建矩阵利用冒号表达式可以线性等间距地建立一个向量来创建矩阵

一般格式是：e1:e2:e3

其中e1为初始值，e2为步长，e3为终止值。

或者为：（start:step:end）例》a=[1:2:10]a=13579第16页，共81页，2022年，5月20日，23点23分，星期四5.矩阵的修改直接修改可用键找到所要修改的矩阵，用键移动到要修改的矩阵元素上即可修改。指令修改可以用A(,)=来修改。第17页，共81页，2022年，5月20日，23点23分，星期四例如a=[120;305;789]a=120305789a(3,3)=0a=120305780第18页，共81页，2022年，5月20日，23点23分，星期四把Matlab工作空间中一些有用的数据长久保存下来的方法是生成mat数据文件。

save——将工作空间中所有的变量存到matlab.mat文件中。二、数据的保存与获取默认文件名第19页，共81页，2022年，5月20日，23点23分，星期四savedata——将工作空间中所有的变量存到data.mat文件中。savedataab——将工作空间中a和b变量存到data.mat文件中。

下次运行Matlab时即可用load指令调用已生成的mat文件。第20页，共81页，2022年，5月20日，23点23分，星期四load——loaddata——loaddataab——mat文件是标准的二进制文件，还可以ASCII码形式保存。即可恢复保存过的所有变量第21页，共81页，2022年，5月20日，23点23分，星期四矩阵加、减（＋,－）运算规则：相加、减的两矩阵必须有相同的行和列两矩阵对应元素相加减。允许参与运算的两矩阵之一是标量。标量与矩阵的所有元素分别进行加减操作。三、矩阵运算第22页，共81页，2022年，5月20日，23点23分，星期四2.矩阵乘（）运算规则：A矩阵的列数必须等于B矩阵的行数标量可与任何矩阵相乘a=[123;456;780];b=[1;2;3];c=a*bc=143223

第23页，共81页，2022年，5月20日，23点23分，星期四d=[-1;0;2];f=pi*df=-3.141606.2832

矩阵除的运算在线性代数中没有，有矩阵逆的运算，在matlab中有两种矩阵除运算。第24页，共81页，2022年，5月20日，23点23分，星期四两种除法：\和/，分别表示左除和右除。如果A矩阵是非奇异方阵，则A\B和B/A运算可以实现。A\B等效于A的逆左乘B矩阵，而B/A等效于A矩阵的逆右乘B矩阵。对于矩阵来说，左除和右除表示两种不同的除数矩阵和被除数矩阵的关系。对于矩阵运算，一般A\B≠B/A。第25页，共81页，2022年，5月20日，23点23分，星期四a^p——a自乘p次幂

方阵>1的整数3.矩阵乘方——a^n,a^p,p^a对于p的其它值,计算将涉及特征值和特征向量，如果p是矩阵，a是标量a^p使用特征值和特征向量自乘到p次幂；如a,p都是矩阵，a^p则无意义。第26页，共81页，2022年，5月20日，23点23分，星期四a=[1,2,3;4,5,6;7,8,9];a^2ans=303642668196102126150※当一个方阵有复数特征值或负实特征值时，非整数幂是复数阵。第27页，共81页，2022年，5月20日，23点23分，星期四a^0.5ans=

0.4498+0.7623i0.5526+0.2068i0.6555-0.3487i1.0185+0.0842i1.2515+0.0228i1.4844-0.0385i1.5873-0.5940i1.9503-0.1611i2.3134+0.2717i第28页，共81页，2022年，5月20日，23点23分，星期四inv——矩阵求逆det——行列式的值eig——矩阵的特征值diag——对角矩阵’——矩阵转置sqrt——矩阵开方4.矩阵的其它运算第29页，共81页，2022年，5月20日，23点23分，星期四5.矩阵的范数矩阵范数的函数为：

(1)norm(V)或norm(V,2)：计算矩阵V的

2—范数。

(2)norm(V,1)：计算矩阵V的1—范数。

(3)norm(V,inf)：计算矩阵V的∞—范数。第30页，共81页，2022年，5月20日，23点23分，星期四6.矩阵的一些特殊操作矩阵的变维

a=[1:12];b=reshape(a,3,4)c=zeros(3,4);c(:)=a(:)矩阵的变向

rot90:旋转;fliplr:左右翻;flipud:上下翻矩阵的抽取

diag:抽取主对角线;(对于非方阵的情况?)

tril:抽取主下三角;

triu:抽取主上三角，然后其余补零元素矩阵的扩展第31页，共81页，2022年，5月20日，23点23分，星期四关系运算

关系符号意义<<=>>===~=小于小于或等于大于大于或等于等于不等于第32页，共81页，2022年，5月20日，23点23分，星期四

关系运算符的运算法则当两个比较量是标量时，直接比较两数的大小。若关系成立，关系表达式结果为1，否则为0。当参与比较的量是两个维数相同的矩阵时，比较是对两矩阵相同位置的元素按标量关系运算规则逐个进行，并给出元素比较结果。最终的关系运算的结果是一个维数与原矩阵相同的矩阵，它的元素由0或1组成。第33页，共81页，2022年，5月20日，23点23分，星期四(3)当参与比较的一个是标量，而另一个是矩阵时，则把标量与矩阵的每一个元素按标量关系运算规则逐个比较，并给出元素比较结果。最终的关系运算的结果是一个维数与原矩阵相同的矩阵，它的元素由0或1组成。注意：其书写方法与数学中的不等式符号不尽相同。第34页，共81页，2022年，5月20日，23点23分，星期四

数组运算指元素对元素的算术运算，与通常意义上的由符号表示的线性代数矩阵运算不同。数组加减(.+,.-)a.+ba.-b7.矩阵的数组运算对应元素相加减（与矩阵加减等效）第35页，共81页，2022年，5月20日，23点23分，星期四2.数组乘除(，./，.\）ab——a，b两数组必须有相同的行和列两数组相应元素相乘。a=[123;456;789];b=[246;135;7910];a.*bans=281841530497290第36页，共81页，2022年，5月20日，23点23分，星期四a=[123;456;789];b=[246;135;7910];a*bans=253746558510985133172第37页，共81页，2022年，5月20日，23点23分，星期四a./b=b.\aa.\b=b./aa./b=b.\a—都是a的元素被b的对应元素除，“/”是斜杠a.\b=b./a—都是b的元素被a的对应元素除，“\”是反斜杠例:a=[123];b=[456];c1=a.\b;c2=b./ac1=4.00002.50002.0000c2=4.00002.50002.0000——给出a,b对应元素间的商.第38页，共81页，2022年，5月20日，23点23分，星期四3.数组乘方(.^)—元素对元素的幂例:a=[123];b=[456];z=a.^2z=1.004.009.00z=a.^bz=1.0032.00729.00(1.^42.^53.^6)第39页，共81页，2022年，5月20日，23点23分，星期四matlab语言把多项式表达成一个行向量，该向量中的元素是按多项式降幂排列的。

f(x)=a0xn+a1xn-1+……an-1x+an

可用行向量p=[a0a1……an-1an]表示poly——产生特征多项式系数向量特征多项式一定是n+1维的特征多项式第一个元素一定是1四、多项式运算第40页，共81页，2022年，5月20日，23点23分，星期四例:a=[123;456;780];p=poly(a)p=1.00-6.00-72.00-27.00p是多项式p(x)=x3-6x2-72x-27的系数matlab描述方法，我们可用：p1=poly2str(p,‘x’)—函数文件，显示数学多项式的形式p1=x^3-6x^2-72x–27注意：多项式中缺少的幂次用‘0’补齐。第41页，共81页，2022年，5月20日，23点23分，星期四2.roots——求多项式的根a=[123;456;780];p=poly(a)p=1.00-6.00-72.00-27.00r=roots(p)---------求由p构成的多项式的根r=12.12-5.73——显然r是矩阵a的特征值

-0.39第42页，共81页，2022年，5月20日，23点23分，星期四当然我们可用poly令其返回多项式形式(这是poly的第二个功能)p2=poly(r)p2=1.00-6.00-72.00-27.00matlab规定多项式系数向量用行向量表示，一组根用列向量表示。第43页，共81页，2022年，5月20日，23点23分，星期四P=poly(r),输入r是多项式所有根，返回值为代表多项式的行向量形式。P=poly(A)，输入是N*N的方阵，返回值p是长度为N+1的行向量多项式，它是矩阵A的特征多项式，也就是说多项式p的根是矩阵A的特征值。第44页，共81页，2022年，5月20日，23点23分，星期四求根的另一种方法str1='x^3-6x^2-72x-27';p1=str2poly(str1);r=roots(p1);注：str2poly实现把一个字符串表示的多项式转换为一个行向量表示的多项式。

poly2str同理。第45页，共81页，2022年，5月20日，23点23分，星期四3.conv多项式乘运算(向量卷积)例:a(x)=x2+2x+3;b(x)=4x2+5x+6;c=(x2+2x+3)(4x2+5x+6)a=[123];b=[456];c=conv(a,b)或c=conv([123],[456])c=4.0013.0028.0027.0018.00p=poly2str(c,‘x’)其中x表示自变量p=4x^4+13x^3+28x^2+27x+18第46页，共81页，2022年，5月20日，23点23分，星期四4.deconv多项式除运算(解卷积)a=[123];c=[4.0013.0028.0027.0018.00]d=deconv(c,a)d=4.005.006.00[d,r]=deconv(c,a)余数c除a后的整数它们之间的关系为:c=conv(a,d)+r第47页，共81页，2022年，5月20日，23点23分，星期四5.多项式导数或微分matlab提供polyder函数计算多项式的导数。命令格式：polyder(p):求p的导数polyder(a,b):求多项式a,b乘积的导数[p,q]=polyder(a,b):求多项式a除以b的商的导数，并以p/q的格式表示。

第48页，共81页，2022年，5月20日，23点23分，星期四例：a=[12345];poly2str(a,'x')ans=x^4+2x^3+3x^2+4x+5b=polyder(a)b=4664poly2str(b,'x')ans=4x^3+6x^2+6x+4第49页，共81页，2022年，5月20日，23点23分，星期四6.多项式的积分matlab提供polyint函数计算多项式的积分。命令格式：polyint(p,k):求多项式p的积分，设积分的常数项为k，polyint(p)默认k=0例：a=[12345];poly2str(a,'x')ans=x^4+2x^3+3x^2+4x+5b=polyint(a,8)b=0.20.51.02.05.08.0poly2str(b,'x')ans=0.2x^5+0.5x^4+x^3+2x^2+5x+8第50页，共81页，2022年，5月20日，23点23分，星期四五、代数方程组求解matlab中有两种除运算左除和右除。对于方程ax=b，a为an×m矩阵，有三种情况：当n=m时，此方程成为“恰定”方程当n>m时，此方程成为“超定”方程当n<m时，此方程成为“欠定”方程

matlab定义的除运算可以很方便地解上述三种方程第51页，共81页，2022年，5月20日，23点23分，星期四1.恰定方程组的解方程ax=b(a为非奇异)x=a-1

b

矩阵逆两种解:x=inv(a)b—采用求逆运算解方程x=a\b—采用左除运算解方程注：若a为奇异的，则Matlab适当给出警告信息或者给出结果为inf。第52页，共81页，2022年，5月20日，23点23分，星期四方程ax=ba=[12;23];b=[8;13];x=inv(a)*b

x=a\bx=x=2.002.003.003.00

=

ax=b例:x1+2x2=82x1+3x2=13第53页，共81页，2022年，5月20日，23点23分，星期四2.超定方程组的解方程的个数大于未知量个数时，方程一般无解。方程解(a'a)x=a'bx=(a‘

a)-1a’b——求逆法(也用到了最小二乘解的原理)x=a\b——matlab用最小二乘法找一个准确地基本解。第54页，共81页，2022年，5月20日，23点23分，星期四例:x1+2x2=12x1+3x2=23x1+4x2=3a=[12;23;34];b=[1;2;3];

解1x=a\b

解2x=inv(a'a)a'b

x=x=1.001.0000.00

=

ax=b第55页，共81页，2022年，5月20日，23点23分，星期四3.欠定方程组的解

当方程数少于未知量个数时,即不定情况,有无穷多个解存在。matlab可求出两个解：用除法求的解x是具有最多零元素的解是具有最小长度或范数的解，这个解是基于伪逆pinv求得的。第56页，共81页，2022年，5月20日，23点23分，星期四x1+2x2+3x3=12x1+3x2+4x3=2a=[123;234];b=[1;2];x=a\bx=pinv(a)bx=x=1.000.8300.330-0.17=ax=b第57页，共81页，2022年，5月20日，23点23分，星期四六、微分方程求解微分方程求解的仿真算法有多种，常用的有Euler（欧拉法）、RungeKutta(龙格-库塔法。Euler法称一步法，用于一阶微分方程第58页，共81页，2022年，5月20日，23点23分，星期四当给定仿真步长时：所以

yn+1=yn+h·f(xn,yn)n=0,1,2…y(x0)=y0第59页，共81页，2022年，5月20日，23点23分，星期四RungeKutta法龙格-库塔法：实际上取两点斜率的平均斜率来计算的，其精度高于欧拉算法。龙格-库塔法：ode23ode45

k1=hf(xn,yn)k2=hf(xn+h,yn+k1)第60页，共81页，2022年，5月20日，23点23分，星期四例：x+(x2-1)x+x=0为方便令x1=x，x2=x分别对x1,x2求一阶导数，整理后写成一阶微分方程组形式

x1=x2x2=x2(1-x12)-x1建立m文件解微分方程······第61页，共81页，2022年，5月20日，23点23分，星期四建立m文件functionxdot=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=[00.25]';[t,x]=ode23('wf',t0,tf,x0)plot(t,x),figure(2),plot(x(:,1),x(:,2))第62页，共81页，2022年，5月20日，23点23分，星期四命令格式:[T,Y]=ODE23(ODEFUN,TSPAN,Y0)建立m文件functiondxdt=wf(t,x)dxdt=[x(2);x(2)*(1-x(1)^2)-x(1)];求解微分方程[t,x]=ode23(@wf,[030],[00.25]);plot(t,x);figure(2)plot(x(:,1),x(:,2))第63页，共81页，2022年，5月20日，23点23分，星期四第64页，共81页，2022年，5月20日，23点23分，星期四七、函数优化寻优函数：fmin——单变量函数fmins——多变量函数constr——有约束条件无约束条件第65页，共81页，2022年，5月20日，23点23分，星期四例1：f(x)=‘x2+3x+2’在[-55]区间的最小值f=fmin('x^2+3*x+2',-5,5)例2：f(x)=100(x2-x12)2+(a-x1)2在x1=a,x2=a2处有最小值functionf=xun(x,a)f=100*(x(2)-x(1).^2).^2+(a-x(1)).^2;x=fmins('xun',[0,0],[],[],sqrt(2))第66页，共81页，2022年，5月20日，23点23分，星期四八、数据分析数据分析相关的函数位于目录：toolboxs\matlab\datafun下Matlab对矩阵操作的规定：如果是向量，则对数据整体操作；如果是矩阵，则对矩阵的列操作。max——各列最大值mean——各列平均值sum——各列求和std——各列标准差var——各列方差sort——各列递增排序cumsum——元素累计和cumprod——元素累计积第67页，共81页，2022年，5月20日，23点23分，星期四八、数据分析数据分析相关的函数位于目录：toolboxs\matlab\datafun下。Matlab对矩阵操作的规定：如果是向量，则对数据整体操作；如果是矩阵，则对矩阵的列操作。max——各列最大值mean——各列平均值sum——各列求和std——各列标准差var——各列方差sort——各列递增排序cumsum——元素累计和cumprod——元素累计积第68页，共81页，2022年，5月20日，23点23分，星期四例：x=[1324],y=[1238;5674]sort(x),sort(y),max(y)>>1234>>12345678>>5678第69页，共81页，2022年，5月20日，23点23分，星期四九、拟合与插值多项式拟合采用最小二乘法对给定的数据进行多项式拟合，最后给出多项式的系数。

p=polyfit(x,y,n)，采用n次多项式p来拟合数据x和y，从而使得y与p(x)最小均方差最小。第70页，共81页，2022年，5月20日，23点23分，星期四x0=0:0.1:1;y0=[-.4471.9783.115.255.024.664.014.583.455.359.22];p=polyfit(x0,y0,3)p=56.6915-87.117440.0070-0.9043xx=0:0.01:1;yy=polyval(p,xx);plot(xx,yy,'-b',x0,y0,'or')第71页，共81页，2022年，5月20日，23点23分，星期四曲线拟合图形用户接口Matlab7.0提供了支持曲线拟合的图形用户接口。在“Figure”窗口“Tools\BasicFitting”菜单中。为了使用该工具，先用待拟合的数据画图。x=0:0.2:10;y=0.25*x+20*sin(x);plot(x,y,'ro');在复选框“Plotfits”中选择“cubic”。第72页，共81页，2022年，5月20日，23点23分，星期四2.插值插值的定义——是对某些集合给定的数据点之间函数的估值方法。当不能很快地求出所需中间点的函数时，插值是一个非常有价值的工具，它可以在已知数据之间寻找估计值，常用