MATLAB与数学实验 第3版 课件汇 07 常微分方程的求解- 13 因子分析法

常微分方程的求解符号解法符号解法数值解法应用常微分方程求解设微分方程初值问题：命令形式1：dsolve(’equation’,’var’)

命令形式2：dsolve(’equation’,’cond1,cond2,…’,’var’)

分为：符号解法和数值解法

符号解法数值解法应用常微分方程求解dsolve('Dy=y^2','x')自变量大写ans=0-1/(C2+x)>>dsolve('Dy=y-1','x')

ans=

C8*exp(x)+1符号解法数值解法应用常微分方程求解求解dsolve('x*D2y-3*Dy=x^2','y(1)=0,y(5)=0','x')大写对应求导阶数conditionans=-1/3*x^3+125/468+31/468*x^4符号解法数值解法应用常微分方程求解数值解法“常微分方程初值问题数值解”的提法：假设(1)式的解存在且唯一……..(1)

不求解析解（无解析解或求解困难）而在一系列离散点求的近似值通常取等步长hxn=x0+nh实验原理

符号解法数值解法应用常微分方程求解数值解法欧拉方法梯形公式……..(1)

对(1)在区间[xi,xi+1]上积分左矩形公式向前欧拉公式……..(*)

欧拉方法用不同的方法求（*）右边定积分各种数值方法实验原理

符号解法数值解法应用常微分方程求解数值解法欧拉方法梯形公式……..(*)

梯形公式梯形公式实验原理符号解法数值解法应用常微分方程求解数值解法例7-44：考虑初值问题：求数值解，并与精确解做比较。精确解为：functionyp=funst(t,y)yp=sec(t)+y*tan(t);解：1)编写函数文件funst.mh=0.05;t=0:h:1;n=length(t);y1=zeros(1,n);y1(1)=pi/2;fori=1:n-1y1(i+1)=y1(i)+h*funst(t(i),y1(i));endyy=(t+pi/2)./cos(t);[t’,y1’,yy’]plot(t,y1,t,yy,'o')2)主程序：实验程序向前欧拉公式符号解法数值解法应用常微分方程求解数值解法实验程序0

1.57081.57080.05001.62081.62280.10001.67491.67920.15001.73361.74030.20001.79721.80680.25001.86651.87920.30001.94191.95830.35002.02432.04480.40002.11442.13970.45002.21342.24420.50002.32242.35970.55002.44282.48770.60002.57642.63020.65002.72512.78970.70002.89152.96900.75003.07863.17180.80003.29033.40290.85003.53153.66800.90003.80833.97480.95004.12874.33361.00004.50334.7581数值解

精确解符号解法数值解法应用常微分方程求解数值解法例7-44：考虑初值问题：试求数值解，并与精确解做比较。精确解为：functionyp=funst(t,y)yp=sec(t)+y*tan(t);解：1)编写函数文件funst.mt0=0;tf=1;y0=pi/2;[t,y]=ode23('funst',[t0,tf],y0);yy=(t+pi/2)./cos(t);plot(t,y,'-',t,yy,'o')[t,y,yy]2)主程序：实验程序龙格库塔方法符号解法数值解法应用常微分方程求解数值解法ans=01.57081.57080.10001.67921.67920.20001.80681.80680.30001.95831.95830.40002.13972.13970.50002.35962.35970.60002.63012.63020.70002.96892.96900.80003.40273.40290.90003.97453.97481.00004.75734.7581t

数值解(y)

精确解(yy)数值解和精确解实验程序符号解法数值解法应用常微分方程求解数值解法实验程序例7-45用数值积分的方法求解微分方程：设初始时间；终止时间初始条件分析：

令原微分方程化为：则：符号解法数值解法应用常微分方程求解数值解法实验程序写成矩阵形式为

functionxdot=exf(t,x)u=1-(t.^2)/(2*pi);xdot=[0,1;-1,0]*x+[0;1]*u;1)编写函数文件exf.m符号解法数值解法应用常微分方程求解数值解法实验程序functionxdot=exf(t,x)u=1-(t.^2)/(2*pi);xdot=[0,1;-1,0]*x+[0;1]*u;1)编写函数文件exf.mt0=0;tf=3*pi;x0=[0;0];2)主程序如下：初始和终止时间初始条件[t,x]=ode23('exf',[t0,tf],x0);y=x(:,1);%用ode23求出x后，取第一列，即为数值解y2=-1/2*(-2*pi-2+t.^2)/pi-(pi+1)/pi*cos(t);%精确解plot(t,y,'b-',t,y2,'go')解析解为：dsolve('D2y+y=1-t^2/(2*pi)','y(0)=0,Dy(0)=0','t')ans=-1/2*(-2*pi-2+t^2)/pi-(pi+1)/pi*cos(t)dy=x(:,2);%y的一阶导数holdon;plot(t,dy,'m--')符号解法数值解法应用常微分方程求解数值解法实验程序例7-46求描述振荡器的VanderPol方程：解：

令则一阶常微分方程为：(状态方程)

符号解法数值解法应用常微分方程求解数值解法实验程序解：一阶常微分方程为：(状态方程)functionxprime=verderpol(t,x)globalmu;xprime=[x(2);mu*(1-x(1)^2)*x(2)-x(1)];1)编写函数文件verderpol.mglobalmu;mu=4;x0=[1;0];[t,x]=ode45('verderpol',[0,20],x0);subplot(1,2,1);plot(t,x);subplot(1,2,2);plot(x(:,1),x(:,2));2)主程序：符号解法数值解法应用常微分方程求解数值解法实验程序例7-47某非刚性物体的运动方程为：取试绘制系统相平面图。初始条件分别对初始值、各个参数的取值做微小的扰动，观察结果符号解法数值解法应用作业应用

符号解法数值解法应用作业应用5、画出卫星的运动轨迹图像。4、完成例7-47线性代数及其应用目录CONTENTS线性代数基本运算案例2：图像压缩案例1：动物繁殖规律k是一个数，A是一个矩阵k*AA\BAX=B,X=A-1B,A必须是方阵数乘矩阵的左除矩阵的右除A/BXB=A,X=AB-1,B必须是方阵矩阵的行列式det(A)A必须为方阵矩阵的逆inv(A)A必须为方阵，矩阵的乘幂A^nA必须为方阵，n是正整数矩阵行变换化简rref(A)求A阶梯形的行最简形式矩阵的基本运算矩阵的转置A’矩阵的秩rank(A)矩阵的特征值、特征向量、特征多项式[V,D]=eig(A)p=poly(A)若A为矩阵，则p为A的特征多项式系数；poly2str(p,’x’)得到多项式的习惯形式A=[1,-1;2,4];[V,D]=eig(A)V=-985/13931292/2889985/1393-2584/2889方阵A的特征向量矩阵D=003方阵A的特征值矩阵A=[1,-1;2,4];p=poly(A)poly2str(p,’x’)p=[1-56]x^2-5x+6在线性代数中用消元法求线性方程组的通解的过程为：1、用初等变换化线性方程组为阶梯形方程组，把最后的恒等式“0=0”去掉；2、如果剩下的方程当中最后的一个等式是等于非零的数，那么方程无解。否则有解；3、在有解的情况下：如果阶梯形方程组中方程的个数r等于未知量的个数，那么方程组有唯一的解；如果阶梯形方程组中方程的个数r小于未知量的个数，那么方程组有无穷多个解。初等变换法解线性方程组例:求解齐次线性方程组的通解Matlab命令为：A=[1-8102;245-1;386-2];rref(A)ans=1.000004.0000001.0000-0.7500-0.25000000分析：∵矩阵A的秩=2＜未知量的个数=4∴方程组有无穷多解同解方程组：取基础解系：通解：为任意实数初等变换法解线性方程组目录CONTENTS线性代数基本运算案例2：图像压缩案例1：动物繁殖规律例:某农场饲养的动物所能达到的最大年龄为15岁，将其分为三个年龄组：第一组：0~5岁；第二组：6~10岁；第三组：11~15岁。动物从第二个年龄组起开始繁殖后代，经长期统计：第二个年龄组的动物在其年龄段平均繁殖4个后代，第三个年龄组的动物在其年龄段平均繁殖3个后代，第一年龄组和第二年龄组的动物能顺利进入下一个年龄组的存活率分别为1/2和1/4。假设农场现有三个年龄段的动物各1000头，问5年，10年及15年后农场饲养的动物总数及农场三个年龄段的动物各将达到多少头？0-5岁6-10岁11-15岁初始状态=1000=1000=1000第一个周期(5年后)

第二个周期(10年后)第三个周期(15年后)案例1：动物繁殖的规律分析：令为0-5岁的动物数,为6-10岁的动物数,

为第i个年龄组在第k个周期的数目。（k=1，2，3）为11-15岁动物数，（k=1，2，3）写成矩阵的形式，有如下矩阵递推关系式：传递矩阵案例1：动物繁殖的规律写成矩阵的形式，有如下矩阵递推关系式：x0=[1000;1000;1000];%初始各年龄组的动物数L=[043;1/200;01/40];x1=L*x0%5年后各年龄组的动物数x2=(L^2)*x0%10年后各年龄组的动物数x3=(L^3)*x0%15年后各年龄组的动物数subplot(1,3,1),pie(x1),title('第一周期后')subplot(1,3,2),pie(x2),title('第二周期后')subplot(1,3,3),pie(x3),title('第三周期后')legend('0-5岁','6-10岁','11-15岁')运行M-文件，得结果x1=%5年后各年龄组的动物数

7000500250x2=%10年后各年龄组的动物数

27503500125x3=%15年后各年龄组的动物数

143751375875案例1：动物繁殖的规律运行M-文件，得结果x1=%5年后各年龄组的动物数

7000500250x2=%10年后各年龄组的动物数

27503500125x3=%15年后各年龄组的动物数

143751375875案例1：动物繁殖的规律

有两家公司R和S经营同类的产品，它们相互竞争。每年R公司保有1/4的顾客，而3/4转移向S公司；每年S公司保有2/3的顾客，而1/3转移向R公司。当产品开始制造时R公司占有3/5的市场份额，而S公司占有2/5的市场份额。问两年后，两家公司所占的市场份额变化怎样，五年以后会怎样？十年以后如何？思考：市场份额如何能够稳定呢？是否所有市场初始分配份额，在经过若干年后均会趋于稳定状态?商品的市场占有率问题商品市场占有率问题

商品的市场占有率问题

clcA=[1/4,1/3;3/4,2/3];x0=[3/5;2/5];fork=1:10disp(['第',num2str(k),'年后两家公司的市场份额为：'])xk=A^k*x0endxi=0.30770.6923[v,d]=eig(A)v=-0.7071-0.40610.7071-0.9138d=-0.0833001.0000思考：最大特征值1对应的特征向量与达到稳态后的xi有什么关系呢？商品的市场占有率问题如何研究长期趋势？

设，为对应于的特征向量，则线性无关。

因此可由线性表示，即存在使得

所以：

因此，当时，

结论：当A的最大特征值为1时，两家公司的长期趋势为最大特征值对应特征向量的倍数。马尔可夫模型概率转移矩阵1、最大特征值为12、其余特征值绝对值小于1一定能达到稳态目录CONTENTS线性代数基本运算案例2：图像压缩案例1：动物繁殖规律案例2：图像压缩问题引入：同学们都上传过照片，是不是经常遇到要求上传照片不能大于100K？如何将手机存储的图像进行压缩，以节省存储空间呢？图像的数据表示。如何将一副彩色图像转化为可以进行运算处理的数字形式。灰度图的处理。灰度图其实就是常见的黑白图片，颜色通道只用一维即可。彩色图像的压缩处理

第一步，通道分离。

第二步，矩阵压缩。

第三步，图像重建。图像的数据表示

一个二维图像，在计算机中通常为一个二维数组，或者一个M*N的二维矩阵（其中，M,N分别为图像的行数和列数）

案例2：图像压缩灰度图像表示

灰度图像也称单色图像，通常由一个二维数组表示一副图像，8位表示一个像素，0表示黑色，255表示白色，1~254表示不同的深浅灰色。img矩阵案例2：图像压缩彩色图像表示RGB图像实一种彩色图形的表示方法，利用3个大小相同的二维数组表示一个像素，3个数组分别代表R、G、B这3个分量，每个像素的每种颜色分量占8位，由[0,255]中的任意数值表示。IMG矩阵imgBimgRimgG案例2：图像压缩问题简化：压缩一张关于主对角线对称的正方形灰度图片，其对应一个实对称矩阵，我们就可以用我们学习过的实对称矩阵理论（实对称矩阵的对角化）进行压缩了。案例2：图像压缩案例2：图像压缩在上式中越大，在中所占比例越大；反之越小，在中所占的比例越小。那我们能否截取中较大部分的特征值的项来近似呢？

案例2：图像压缩[v,d]=eigs(A,k)

：返回

k

个最大幅值特征值

imread(filenme)：指定的文件读取图像

imshow(I)

在图形中显示图像II

=rgb2gray(RGB)

将真彩色图像

RGB

转换为灰度强度图像

I

案例2：图像压缩clcclearclfIMG=imread('img2.jpg');img=rgb2gray(IMG);forh=1:180

forl=1:himg(l,h)=img(h,l);

end

endsubplot(221),imshow(img)title('原灰度图像')img=double(img)[M,N]=eigs(img,18)imgc1=M*N*M';subplot(222),imshow(uint8(imgc1))title(‘前18个特征值重构图像')[M,N]=eigs(img,36)imgc2=M*N*M';subplot(223),imshow(uint8(imgc2))title(‘前36个特征值重构图像')[M,N]=eigs(img,54)imgc3=M*N*M';subplot(224),imshow(uint8(imgc3))title(‘前54个特征值重构图像')案例2：图像压缩案例2：图像压缩imgr=IMG(:,:,1);imgg=IMG(:,:,2);imgb=IMG(:,:,3);imgc1(:,:,1)=uint8(imgcr1);imgc1(:,:,2)=uint8(imgcg1);imgc1(:,:,3)=uint8(imgcb1);第一步，通道分离第二步，矩阵压缩,同灰度图片

第三步，图像重建案例2：图像压缩问题一般化：奇异值分解

案例2：图像压缩奇异值分解

案例2：图像压缩提取前k项作为A的近似，则：

m行n列

奇异值分解

案例2：图像压缩那么如何求出SVD分解后的U,Σ,V这三个矩阵呢？

将n×n方阵

进行特征分解，得到的特征值和特征向量满足下式：

奇异值分解

将

的所有特征向量张成一个n阶矩阵V，就是SVD公式里的V矩阵。一般将V中的每个特征向量叫做A的右奇异向量。

将m×m方阵进行特征分解，得到的特征值和特征向量满足下式：

将的所有特征向量张成一个m阶矩阵U，就是SVD公式里的U矩阵。一般将U中的每个特征向量叫做A的左奇异向量。

求出SVD分解后的U,Σ,V这三个矩阵？1、求矩阵V2、求矩阵U

案例2：图像压缩奇异值分解

3、求奇异值矩阵案例2：图像压缩Σ2：矩阵或的特征值组成的对角矩阵Σ：矩阵A的奇异值组成的对角矩阵因此，矩阵或矩阵特征值和矩阵A的奇异值满足：[U,S,V]

=svds(A,k)

返回

k

个最大奇异值的左奇异矢量U、奇异值的对角矩阵S以及右奇异量V

案例2：图像压缩案例2：图像压缩作业1、已知矩阵A=,实现下列操作：（1）添加零元素，使之成为3×3的方阵；（2）在以上操作的基础上，将第三行元素换为（135）；（3）在以上操作的基础上，提取矩阵中第2个元素以及第3行第2列的元素。2、求矩阵A=的特征多项式、特征值和特征向量。3、求线性方程组的通解：4、同学们自行将自己的彩色证件照分别进行灰度压缩，彩色压缩5、课件上商品市场占有率问题概率相关运算目录CONTENTS统计作图常见的概率分布蒙特卡洛模拟常见的概率分布

分布

离散型随机变量连续型随机变量均匀分布二项分布泊松分布几何分布均匀分布指数分布正态分布分布t分布

F分布字符unidbinopoissgeounifexpnormchi2tf功能概率密度分布函数逆概率分布均值与方差随机数生成字符pdfcdfinvstatrnd表8-2概率分布的命令字符表8-3运算功能的命令字符概率密度为y=normpdf(x,mu,sigma)的正态分布的密度函数y=normpdf(x)标准正态分布的密度函数y=normcdf(x,mu,sigma)的正态分布的分布函数y=normcdf(x)标准正态分布的分布函数分布函数为常见的概率分布x=-6:0.01:6;y1=normpdf(x);z1=normcdf(x);y2=normpdf(x,0,2);z2=normcdf(x,0,2);subplot(1,2,1),plot(x,y1,x,y2);legend('N(0,1)','N(0,2^2)')subplot(1,2,2),plot(x,z1,x,z2);legend('N(0,1)','N(0,2^2)')常见的概率分布目录CONTENTS统计作图常见的概率分布蒙特卡洛模拟统计作图输入数据在数据较小、较少的情况下输入Matlab交互环境下输入M文件的形式输入数据数据量较大，且不以计算机可读形式存在load*.M读数据文件的命令读入load*.txt读数据文件的命令读入xlsread(‘*.xlsx’)基本统计函数函数名称功能简介max(x)求最大值min(x)求最小值median(x)求中值range(x)求极差mean(x)求算术平均值std(x)求样本标准差var(x)求样本方差cov(x)求协方差矩阵统计作图例8-19某班（共有120名学生）的高等数学成绩如下：746378768956709789947688658372413972736814764570904654617576495778666474788786734767216679676865568466736872766570945365777853745950986789786392548784806364856669696054753330627465847355857576817183725684767567653594594745677536788294708475根据以上数据作出该门课程成绩的频数表和直方图。统计作图解:(1)数据输入：将以上数据以一列的形式存为A.txt文件,用

(2)用hist命令作频数表和直方图：(区间个数为5，可省略)[N,X]=hist(A,5)120名学生高数成绩的频数表；hist(A,5)120名学生高数成绩的直方图；loadA.txt命令读入数据。频数是如何计算的？统计作图a1=min(A);a2=max(A);disp([‘成绩最小值’,blanks(4),‘最大值'])disp([a1,a2])[N,X]=hist(A,5)N=310226025X=22.400039.200056.000072.800089.6000成绩最小值最大值149884—÷5=16.814+16.8=30.8所以，成绩在14~30.8之间有3人。22.4从哪儿来的？统计作图hist(A,5)%直方图

M=68.958371.500084.0000249.569715.7978例8-20求例8-19中A的均值、中位数、极差、方差和标准差。解：在命令窗口输入：M=[mean(A),median(A),range(A),var(A),std(A)]统计作图目录CONTENTS统计作图常见的概率分布蒙特卡洛模拟蒙特卡洛模拟问题的引入

解决确定性问题计算定积分以及重积分如右图，如何计算不规则图形（阴影区域）的面积？记：事件A——随机点落入阴影区域内总点数点数(阴影)频率概率

大数定律如何生成随机点？横坐标与纵坐标均是随机数蒙特卡洛模拟法：利用随机数进行数值模拟，并获得问题解的方法。

蒙特卡洛模拟假定在[a,b]区间上定义了一个非负连续函数f(x)，则曲边梯形的面积可以表示为：解决确定性问题计算定积分以及重积分，其中对于每个i值，，则这种点数k与总数n之比k/n即表示曲边梯形ABCD与矩形区域总面积之比。随机取n个数对如果点落在曲边梯形ABCD内，实例：计算定积分并与精确值作比较。x=unifrnd(0,pi,1,10000)y=unifrnd(0,pi,1,10000)%x=unifrnd(0,pi,2,100)蒙特卡洛模拟解决确定性问题计算定积分以及重积分蒙特卡洛模拟的优点：精度与维数无关推广计算多重积分unifrnd(a,b,m,n)生成m*n个服从[a,b]均匀分布的随机数实例：计算定积分并与精确值作比较。x=unifrnd(0,pi,1,10000)y=unifrnd(0,pi,1,10000)%x=unifrnd(0,pi,2,100)蒙特卡洛模拟解决确定性问题计算定积分以及重积分unifrnd(a,b,m,n)生成m*n个服从[a,b]均匀分布的随机数clcn=100000;x=unifrnd(0,pi,1,n);y=unifrnd(0,pi,1,n);count=y<sin(x);%判断那些点落在阴影区域内freq=sum(count)/n;%频率f=freq*pi^2%蒙特卡洛估计值ff=quad('sin(x)',0,pi)%复合辛普生估计值试验序号12345672315124222521252418272512492562472512622580.40.60.21.00.20.40.80.440.500.420.480.360.540.5020.4980.5120.4940.5240.5160.500.502抛硬币实验将一枚硬币抛掷5次、50次、500次,各做7遍,观察正面出现的次数及频率.波动最小随n的增大,频率

f呈现出稳定性蒙特卡洛模拟解决随机性问题抛硬币实验将一枚硬币抛掷100次,1000次,10000次

观察正面出现的次数及频率.n=0;m=100;%模拟抛硬币的次数x=unidrnd(2,1,m)-1;y=0;fori=1:m

if(x(i)==1)y=1;

elsey=0;

endn=n+y;endf=n/m蒙特卡洛模拟解决随机性问题unidrnd(N,m,n)功能：生成m*n个最大值为N的服从均匀分布的随机正整数生日问题蒙特卡洛模拟解决随机性问题在100人的团体中，如果不考虑年龄的差异，研究是否有两个以上的人生日相同。假设每人的生日在一年365天中的任意一天是等可能的，那么随机找n个人（n<365）:问这些人生日各不相同的概率是多少？

至少有两个人生日相同的概率是多少？forn=1:100p0(n)=prod(365:-1:365-n+1)/365^n;p1(n)=1-p0(n);endn=1:100;plot(n,p0,n,p1,'--')xlabel(‘人数’),ylabel(‘概率’)legend(‘生日各不相同的概率’,‘至少两人相同的概率’)>>p1(30)ans=0.7063蒙特卡洛模拟解决随机性问题生日模拟实验n=30时，至少有两个人生日相同的概率是多少？0.70631）模拟一个班---------------------------做一次试验一个班30个人---------------------产生30个随机数生日只可能是1~365天中的某一天---------随机数的范围是[1,365]每个人生日在365天的任意一天是等可能的--------服从均匀分布如果随机数中出现相同的数（生日相同），就记为1，否则记为02）模拟100个班--------------------------------------做100次试验最后统计1的个数（出现生日相同的频数）生日相同的概率=频数/1003）与理论概率做比较0.7063蒙特卡洛模拟解决随机性问题生日模拟实验（模拟1个班）蒙特卡洛模拟解决随机性问题生日模拟实验（模拟100个班）作业1、生日模拟问题，计算频率，并画出随着试验次数n的增大，频率和理论概率的关系图2、画出正态分布N(0,0.42),N(0,12),N(-2,22),N(1,22)的概率密度图和分布函数图3、下面的数据是一个专业50名大学新生的测验分数：将这组数据分成6-8组，画出频数直方图，求出样本均值和方差4、设某种药物的临床有效率为0.95，现有10人服用，则至少8人治愈的概率是多少?作业5、模拟掷骰子，计算出现5点的频率，并画出随着试验次数n的增大，频率和概率的关系图或模拟抛硬币，计算正面出现的频率，并画出随着试验次数n的增大，频率和概率的关系图。（unidrnd(N):生成最大值为N的随机正整数。）6、用蒙特卡洛方法计算定积分的近似值，并与精确值

作比较7、用蒙特卡洛方法计算二重积分作业3题数据90,76,69,51,71,40,88,79,68,77,96,69,80,71,86,52,41,60,81,72,92,81,99,77,100,79,66,71,84,73,67,70,86,75,60,80,77,91,93,64,74,76,83,81,83,88,80,92,83,64优化目录CONTENTS简单优化1产销量的最佳安排2再谈拟合3优化01

目标函数标准形式：

优化问题是工程技术、经济管理和科学研究等领域中最常遇到的一类问题

设计师要在满足强度要求等条件下选择材料的尺寸，使结构总质量最轻;公司经理要根据生产成本和市场需求确定产品价格，使所获利润最高;调度人员要在满足物资需求和装载条件下安排从各供应点到各需求点的运量和路线，使运输总费用最低;投资者要选择一些股票、债券"下注"，使收益最大，而风险最小……

优化问题归结为求函数极值问题011、fminbnd函数（求解一元函数优化问题，局部最优）

命令格式：[xmin,ymin]=fminbnd(f,a,b),ymin为极小值[xmin,ymin,exitflag,output]=fminbnd(f,a,b)优化过程的信息：iterations(迭代次数),funcCount(函数计算次数),algorithm(算法：黄金分割法，抛物线法),message(退出消息)返回函数fminbnd的求解状态（成功或失败）xmin=fminbnd(f,a,b),求函数f在区间[a,b]上的极小值，其中f是用来求极值的函数，可以是函数名，也可以是函数表达式。xmin为极小值点。简单优化01

黄金分割法

思考：如何缩小搜索区间的范围？简单优化01

黄金分割法

黄金分割点

简单优化01

黄金分割法的步骤

1.确定初始搜索区间[a,b]

简单优化01例1：求函数的极值。

[xmin,ymin]=fminbnd('x.^3+x.^2-5*x-5',-5,5)运行结果：xmin=1.0000,ymin=-8.0000[x2,y2]=fminbnd('-(x.^3+x.^2-5*x-5)‘,-5,5)ymax=-y2运行结果：极大值点x2=-1.6667,极大值为ymax=1.4815因此，函数y在[-5,5]上的极小值为-8，极大值为1.4815。

简单优化简单优化01[xmin,ymin]=fminbnd('sin(2*x+1)+2*sin(4*x+3)',-3,3)运行结果：xmin=1.9832,ymin=-2.9640求所有极小值点：例2：求函数在区间的极小值。

[xmin1,ymin1]=fminbnd('sin(2*x+1)+2*sin(4*x+3)',-1.5,-1)[xmin2,ymin2]=fminbnd('sin(2*x+1)+2*sin(4*x+3)',0,1)[xmin3,ymin3]=fminbnd('sin(2*x+1)+2*sin(4*x+3)',1.5,2.3)[xmin4,ymin4]=fminbnd('sin(2*x+1)+2*sin(4*x+3)',-3,-2.2)思考：为什么只求出一个极小值？012、fminunc函数（求解一元函数及多元函数的极小值，局部最优）标准形式：

求解的基本思想：迭代法思考：如何快速找到最低点(谷底)？1、初值简单优化012、fminunc函数标准形式：

2、方向求解的基本思想：迭代法思考：如何快速找到最低点(谷底)？1、初值简单优化012、fminunc函数标准形式：

3、步长第k步的搜索方向第k步的步长方向不同，或步长不同，对应着不同的数值方法。求解的基本思想：迭代法2、方向思考：如何快速找到最低点(谷底)？1、初值常见的方法有：梯度下降法、最速下降法、牛顿法、拟牛顿法等。简单优化012、fminunc函数梯度下降法：负梯度方向初始点：初始点：简单优化012、fminunc函数梯度下降法：负梯度方向步长：……简单优化012、fminunc函数梯度下降法：负梯度方向03611.83.621.082.1630.6481.29640.38880.777650.233280.4665660.1399680.27993670.08398080.167961680.050388480.1007769690.0302330880.060466176100.0181398530.036279706110.0108839120.021767823

简单优化01……简单优化01……简单优化012、fminunc函数[x,fval]=fminunc(fun,X0)[x,fval,exitflag,output]=fminunc(fun,X0)

functionf=fun(x)f=2*x(1)^2+x(2)^2;[x,fval,exitflag,output]=fminunc('fun',[3,3])

简单优化01例4：求函数的极小值。functionf=fun1(x)f=2*x(1)^2-x(1)^4+x(1)^6/6-x(1)*x(2)+x(2)^2;[x_1,minf_1]=fminunc('fun1',[-1.5,-1])[x_2,minf_2]=fminunc('fun1',[0,-0.5])[x_3,minf_3]=fminunc('fun1',[1.5,0.5])x_1=-1.6453-0.8227minf_1=0.7155x_2=1.0e-07*

0.1072-0.0466minf_2=3.0137e-16x_3=1.64530.8227minf_3=0.7155简单优化目录CONTENTS无约束优化1产销量的最佳安排2再谈拟合3产销量的最佳安排02某厂生产一种产品有甲、乙两个牌号，讨论在产销平衡的情况下如何确定各自的产量，使总利润最大.所谓产销平衡指工厂的产量等于市场上的销量。问题分析总利润既取决于销量和价格，也依赖于产量和成本。

甲的价格甲的成本甲的销量乙的价格乙的成本乙的销量产销量的最佳安排02模型假设1、价格与销量成线性关系影响甲的价格p1的因素：1）随甲的销量x1的增长而降低2）乙的销量x2的增长也会使甲的价格有稍微的下降可以简单地假设价格与销量成线性关系，即：

同理：

总利润既取决于销量和价格，也依赖于产量和成本。产销量的最佳安排02模型假设1、价格与销量成线性关系

2、成本与产量成负指数关系

甲的成本q1随其产量(等于甲的销量x1)的增长而降低，且有一个渐进值，可以假设为负指数关系，即：同理：

产销量的最佳安排02模型建立

总利润：若根据大量的统计数据，求出系数b1=100,a11=1,a12=0.1,b2=280,a21=0.2,a22=2,r1=30,λ1=0.015,c1=20,r2=100,λ2=0.02,c2=30

则问题转化为无约束优化问题：求甲，乙两个牌号的产量x1,x2，使总利润z最大。产销量的最佳安排02模型求解

用迭代法求解，先估计初始值：简化模型，先忽略成本，并令,a12=0,a21=0，则问题转化为求

的极值显然其解为x1=b1/2a11=50，x2=b2/2a22=70，我们把它作为原问题的初始值。产销量的最佳安排02编程求解1.建立函数文件：functionz_minus=fun(x)y1=((100-x(1)-0.1*x(2))-(30*exp(-0.015*x(1))+20))*x(1);y2=((280-0.2*x(1)-2*x(2))-(100*exp(-0.02*x(2))+30))*x(2);z_minus=-y1-y2;2.主程序：x0=[50,70];xmin=fminunc(‘fun’,x0)z_minus=fun(xmin)

3.计算结果：x=23.9025,62.4977,z_minus=-6.4135e+003即甲的产量为23.9025，乙的产量为62.4977，最大利润为6413.5。目录CONTENTS无约束优化1产销量的最佳安排2再谈拟合3再谈拟合03人口模型——Malthus模型年份人口统计数字年份人口统计数字19718.5229198110.007219728.7177198210.165419738.9211198310.300819749.0859198410.435719759.2420198510.585119769.3717198610.750719779.4974198710.930019789.6259198811.102619799.7542198911.270419809.8705199011.4333Malthus模型：目标：根据前15年的数据进行拟合，求p0,r，使得

最小即：再谈拟合033、lsqcurvefit函数（least-squarescurve-fitting）X=lsqcurvefit(FUN,X0,XDATA,YDATA)已知的数据点迭代初始点注意：FUN返回FUN(X,XDATA)，而不是sum((FUN(X,XDATA)-YDATA).^2)用于最小二乘法求解非线性曲线拟合问题min

sum{(FUN(X,XDATA)-YDATA).^2}Xfunctionpt=population1(x,tdata)pt=x(1)*exp(x(2)*tdata);%x(1)=p0,x(2)=rFUN(X,XDATA),X是未知系数组成的向量，

Malthus模型中，再谈拟合03tdata=1971:1985;tdata=tdata-1970;pdata=[8.5229,8.7177,8.9211,9.0859,9.2420,9.3717,9.4974,9.6259,9.7542,9.8705,10.0072,10.1654,10.3008,10.4357,10.5851];x0=[1,0.15];%x的初值，迭代初始值

x=lsqcurvefit('population1',x0,tdata,pdata)%极小值点（p0，r）ts=1971:2015;ts=ts-1970;y=population1(x,ts);plot(ts+1970,y)人口模型——Malthus模型再谈拟合03人口模型——阻滞增长模型（Logistic模型）由于自然资源、环境条件等因素对人口的增长起着阻滞作用，并且随着人口的增加，阻滞作用越来越大.1）人口增长率r为当时人口数量p(t)的减函数，最简单的假定：r是固有增长率；2）自然资源和环境条件年容纳的最大人口数量为pm模型假设再谈拟合031）人口增长率r为当时人口数量p(t)的减函数，最简单的假定：r是固有增长率；2）自然资源和环境条件年容纳的最大人口数量为pm当p=pm时，增长率为0，即：所以：代入到Malthus模型模型假设模型建立再谈拟合03用分离变量法得：模型求解两边分别积分，得：所以：functionpt=population2(x,tdata)pt=x(1)./(1+(x(1)/8.5278-1)*exp(-x(2)*tdata));%x(1)=pm,x(2)=r

Logistic模型再谈拟合03模型求解functionpt=population2(x,tdata)pt=x(1)./(1+(x(1)/8.5229-1)*exp(-x(2)*tdata));%x(1)=pm,x(2)=rtdata=1971:1985;tdata=tdata-1970;pdata=[8.5229,8.7177,8.9211,9.0859,9.2420,9.3717,9.4974,9.6259,9.7542,9.8705,10.0072,10.1654,10.3008,10.4357,10.5851];x0=[16,0.15];%x的初值x=lsqcurvefit(‘population2’,x0,tdata,pdata)%最优解（pm，r）ts=1971:2015;ts=ts-1970;y=population2(x,ts);plot(ts+1970,y)

作业043、用下面一组数据拟合

中的参数a，b，k1、分别用命令fminbnd,fminunc求以下函数在区间[-8，8]的所有极小值。2、求peaks函数的最小值>>peaksz=3*(1-x).^2.*exp(-(x.^2)-(y+1).^2)...-10*(x/5-x.^3-y.^5).*exp(-x.^2-y.^2)...-1/3*exp(-(x+1).^2-y.^2)数学建模实例115层次分析方法引例1161、大学毕业生就业选择问题获得大学毕业学位的毕业生,在面临就业的“双向选择”时,用人单位与毕业生都有各自的选择标准和要求就毕业生来说:选择单位的标准和要求是多方面的,例如:能发挥自己的才干为国家做出较好贡献(即工作岗位适合发挥专长);工作收入较好(待遇好）生活环境好(大城市、气候、等工作作条件)单位名声好(声誉)工作环境好(人际关系和谐等)发展晋升机会多(如新单位或单位发展有前途)等生活环境工作选择贡献收入发展声誉工作环境可供选择的单位p1,p2,...,pn117引例2、在基础研究应用研究和数学教育中选择一个领域申报科研课题，要考虑成果的贡献(实用价值、

科学意义),可行性(难度周期和经费)和人才培养，相应的层次结构图如下所示经费科研课题实用价值学术意义人才培养周期基础难度可行性贡献应用教育1183、“五一”假期快到了，张同学决定假期去踏青。他想去桂林、黄山和北戴河三个踏青地点,但由于时

间的原因，他只能在这三个地点中选一个来作为踏青的目的地。相应的层次结构图如下所示选择旅游地景色收费居住饮食旅途可供选择的旅游目的地（桂林、黄山、北戴河）引例119这些问题属于决策问题,即做一件事情有多个选择,怎样才能选择最好的一个。要在多个选择对象中选择其中一个，人们往往要根据自己的目标和有利于目标实现的多种因素的考量来做出最后的选择要达到使决策问题化为具有条理化和层次化的结构模型,需要我们先把复杂问题分解为一些因素,然后把这些因素按其属性及关系形成若干层次，上一层次的因素作为准则对下一层次的有关因素起支配作用，把这些层次可以分为三类：目标层：这一层次中只有一个因素，一般它是决策问题的预定目标或理想结果，处于层次结构的第一层准则层：这一层次中包含了为实现目标所涉及的中间环节，它可以由若干个层次组成，包括所需考虑的准则、子准则，处于层次结构的中间层方案层：这一层次包括了为实现目标可供选择的各种措施、决策方案等，因此也称为措施层，处于层次结构的最底层引例总结120以踏青问题为例建立层次分析方法模型踏青问题可以表示为三层的递阶层次结构。第一层(选择最佳旅游地)是目标层，第二层(判断旅游地的倾向)是准则层，第三层(旅游地点)是方案层，它们之间用线段连接表示它们之间的联系要依据喜好对三个层次相互比较进行综合判断，在三个旅游地中确定哪一个为最佳地点。如下图所示选择旅游地景色收费居住饮食旅途桂林黄山北戴河目标层准则层方案层121决策过程中遇到的困难：在确定影响某因素的诸因子在该因素中所占的比重时，遇到的主要困难是这些比重常常不易定量化。此外，当影响某因素的因子较多时，直接考虑各因子对该因素有多大程度的影响时，常常会因考虑不周全、顾此失彼而使决策者提出与他实际认为的重要性程度不一致的数据，甚至有可能会提出一组隐含矛盾的数据以踏青问题为例建立层次分析方法模型解决方案：建立成对比较矩阵成对比较矩阵122成对比较矩阵定义:假设要比较

个因子对某因素的影响大小，每次取其中的两个因子和以

表示

和

对的影响大小之比，全部比较结果用矩阵表示,为成对比较判断矩阵(简称判断矩阵)态度、喜好等主观因素量化的形式。层次分析法中较多使用1-9标度，通过对两两因素进行比较获得相对重要程度的关系(重要程度也可解释为偏好、可能性等)，这种重要程度用数字1-9及其倒数表示为相互间的倍数。标度概念:

数字1-9标度123心理学家认为成对比较的因素不宜超过9个，即每层不要超过9个因素标度含义1与的影响相同3比影响稍强5比影响强7比

影响明显的强9比

影响绝对的强2,4,6,8

比

影响之比在上述相邻等级之间1,1/2,…,1/9与

影响之比的倒数

踏青问题的成对比较矩阵124项目景色费用饮食居住旅途景色11/2553费用21775饮食1/51/711/21/3居住1/51/7211/2旅途1/31/5321容易看出，若与

对的影响之比为，则与

对

的影响之比为：

如果成对比较矩阵构造准确，则应满足一致性条件

一致性检验125

随机一致性指标一致性指标注意到上述一致性指标

是一个绝对量，不易说明取值多小才算是很小。为确定

的不一致程度的容许范围，还要找出衡量-致性指标的标准，要引人一个相对的量来描述“很小”的取值。为此Saaty借助随机试验的方式引人了随机一致性指标

它描述了一致性指标

的平均值

12345678910000.580.901.121.241.321.411.451.49随机一致性指标

的得出为对每个固定的

，随机的构造100至500个正互反矩阵

，

然后计算每一个矩阵的一致性指标，再取平均值，由此得到随机一致性指标的值如上表所示

126一致性比例

通常0.1被用作临界值若规定

的临界值是0.1，即作为比较矩阵的一致性是可以接受的标准，它暗示了该比较矩阵的一致性指标

比其平均值的10%还小可以认为是一个很小的数了。如果没有通过一致性检验，应该对比较矩阵作适当修正以使其通过一致性检验

一致性检验127求各因素的权重项目景色费用饮食居住旅途景色11/2553费用21775饮食1/51/711/21/3居住1/51/7211/2旅途1/31/5321按列归一化3.73331.98571815.59.8333项目景色费用饮食居住旅途景色0.25840.25180.27780.32260.3051费用0.51680.50360.38890.45160.5085饮食0.05170.07190.05560.03230.0339居住0.05170.07190.11110.06450.0508旅途0.08610.10070.16670.12900.10170.28320.47400.05600.07000.11691.47552.47400.24900.35330.5971

5.0988

0.0247

1.120.0221

踏青问题第2层对第1层一致性检验及权重128项目景色费用饮食居住旅途景色11/2553费用21775饮食1/51/711/21/3居住1/51/7211/2旅途1/31/5321最大特征值：5.0976[V,D]=eigs(A)0.49560.83260.08400.11870.20040.02440.0218

通过一致性检验权重clearclcclfA=[1,1/2,5,5,3;2,1,7,7,5;1/5,1/7,1,1/2,1/3;1/5,1/7,2,1,1/2;1/3,1/5,3,2,1];[V,D]=eigs(A);Eigen_value=real(D)Eigen_value_max=Eigen_value(1)Characteristic_vector=real(V)Characteristic_vector_max=Characteristic_vector(:,1)Normalization_vector_max=Characteristic_vector_max(:,1)./sum(Characteristic_vector_max(:,1))CI=(Eigen_value_max-length(A(:,1)))/(length(A(:,1))-1)CR=CI/1.12踏青问题第2层对第1层一致性检验及权重

标准化后的权重129把如上关系表示为数学形式，有：

类似的，有：

设目标层只有一个因素，第二层包含

个因素,他们关于

的权重分别为；第三层包含

个因素，

；

第三层对第二层的每个因素的权重为

组合权重

组合一致性检验及权重

130组合一致性检验方法：先逐层进行组合一致性检验，然后考虑整个层次逐层进行组合一致性检验

设第

层有

个成对比较矩阵,其对应的

个一致性指标，随机一致性指标为，定义第层的一致性指标

和随机一致性指:

其中,是第层对第

层的组合权向量。用

除以

就得到第层的一致性比率。

组合一致性比率

若整个层的比较判断通过一致性检验

组合一致性检验及权重131景色桂林黄山北戴河桂林125黄山1/212北戴河1/51/21费用桂林黄山北戴河桂林11/31/8黄山311/3北戴河831饮食桂林黄山北戴河桂林113黄山113北戴河1/31/31

踏青问题组合一致性检验及权重

标准权重132居住桂林黄山北戴河桂林134黄山1/311北戴河1/411

旅途桂林黄山北戴河桂林111/4黄山111/4北戴河441

第3层的权重矩阵为：

第3层对第1层组合权向量为：

踏青问题组合一致性检验及权重

标准化后的权重133第3层的一致性指标和随机一致性指标：

第3层对第1层的组合一致性比率：

通过一致性检验结论：通过组合权向量可知方案3(北戴河)在踏青旅游选择占权重为0.458，明显大于方案1(桂林权重为0.294)方案2(黄山权重为0.247)。故在所设定标准的前提下,他们应该去北戴河踏青问题组合一致性检验及权重标准化后的权重标准化后的权重134层次分析方法步骤（1）建立层次结构模型，其中最高层为单个目标层，最底层是要决策的方案；（2）从第2层开始由上到下的顺序对每一层构造出各层次中的所有成对比较矩阵；并计算每一个成对比较矩阵

的最大特征值和特征向量，计算

判别是否成立，若成立，一致性检验通过，将特征向量归一化得到权重，否则，重新构造该成对较矩阵（3）按公式计算组合权重

（4）按公式

计算组合一致性比率（5）按公式

进行组合一致性检验，通过则根据组合权向量的分量做出决策；若不通过，重新考虑模型结构或重新构

造那些一致性比率较大的成对比较矩阵

踏青问题组合一致性检验及权重作业135基于省时、收入、岸间商业、当地商业、建筑就业5项因素拟用层次分析法在建桥梁、修隧道、设渡轮这3个方案中选最佳，画出目标为“越海方案的最优经济效益”的层次结构图，用两种层次分析法选出方案其中，第一种方法交Excel表格即可，第二种方法交Matlab程序文件。主成分分析方法问题的提出

多变量问题是经常会遇到的。变量太多，无疑会增加分析问题的难度与复杂性，而且在许多实际问题中，多个变量之间是具有一定的相关关系的。因此，人们会很自然地想到，能否在相关分析的基础上，用较少的新变量代替原来较多的旧变量，而且使这些较少的新变量尽可能多地保留原来变量所反映的信息？

主成分分析就是设法将原来指标重新组合成一组新的互相无关的几个综合指标来代替原来指标。同时根据实际需要从中可取几个较少的综合指标尽可能多地反映原来的指标的信息。

从数学角度来看，这是一种降维处理技术。目