matlab教学第一章方程(组)模型

1、第一章第一章 方程（组）模型方程（组）模型2022-3-161本章学习目的：本章学习目的：l本章我们就是要学习求解线性方程组、非线性方程（组）的方法，以及利用数学软件和计算机对方程和方程组进行求解。l考虑求方程f(x)=0的解，我们通常采用这样的几种方法：因式分解法、图形放大法、数值迭代逼近法。2022-3-162本章学习目的：本章学习目的：l复习求解方程的基本原理和方法，掌握解方程的图形放大法和迭代算法；l熟练掌握用MATLAB软件的函数来求解方程和方程组；l通过范例展现求解实际问题的初步建模过程和MATLAB程序设计。2022-3-1631.1 图形放大法：l由于计算机的广泛应用，可以非常

2、方便地作出函数函数f(x)的图形的图形（曲线），找出曲线与x轴的交点的横坐标值，就可求出f(x)0的近似根。l这些值尽管不精确，但是直观，方程有多少个根、在什么范围，一目了然。并且可以借助于计算机使用图形局部放大图形局部放大功能，将根定位定位得更加准确准确。2022-3-164用图形放大法求解方程用图形放大法求解方程f(x)0的步骤：的步骤：l建立坐标系，作曲线f(x)；l观察f(x)与x轴的交点；l将其中的一个交点进行局部放大；l该交点的横坐标值就是方程的一个根；l对所有的交点进行相同的处理，就得到方程的所有解。2022-3-165例例1.1 求方程求方程 所有的根及大致分所有的根及大致分布

3、范围，欲寻求其中的一个实根，并且达到一布范围，欲寻求其中的一个实根，并且达到一定的精度。定的精度。l画出 的图形；l程序如下：lx=-6:0.01:6;ly=x.5+2*x.2+4;lplot(x,y)lgrid on; %画坐标格52240 xx42)(25xxxf2022-3-166l我们可以看出方程在22范围有一个实根。-6-4-20246-8000-6000-4000-2000020004000600080002022-3-167逐次缩小范围得到较精确的根。逐次缩小范围得到较精确的根。-2-1.5-1-0.500.511.52-20-1001020304050-1.6-1.58-1.5

4、6-1.54-1.52-1.5-1.48-1.46-1.5-1-0.500.51最后我们可以看出这个实根的值在-1.56-1.54之间。2022-3-1681.2 简单迭代法简单迭代法l迭代算法步骤：迭代算法步骤： l对方程f(x)=0求解l对方程经过简单变形得到 （不是唯一的），x 被称之为不动点；l设置为迭代初值 ，迭代过程为 ，n=0,1,2当两次迭代结果之差小于某个设定的误差值时，我们认为迭代结果是收敛的，可得到结果的近似值 。)(xx)(1nnxx1nxx2022-3-169nx例1.2 求方程求方程 的非负实根。的非负实根。l解：由于函数 连续，并且在x=0和x=1处函数值符号相反

5、，可以判断函数在区间（0，1）必有零点，即方程在（0，1）内必然存在根。l先将函数变形为 ；l设置迭代初值为0，编程进行迭代。03xexxex 33/xex 2022-3-1610ln=1;lx=0;ly=exp(x)/3;lys=vpa(y,10); %ys取取y的的10位有效数字位有效数字lz=abs(y-x);lwhile z10(-5)l x=y;l y=exp(x)/3;l ys=vpa(y,10);l z=abs(y-x);l n=n+1;lendln,y,ys n = 21y = 0.6190ys =.6190471917 从该结果可以看出，迭代从该结果可以看出，迭代21次后次后

6、两次迭代的结果误差值满足小于两次迭代的结果误差值满足小于的条件，结果收敛，迭代结果为的条件，结果收敛，迭代结果为0.6190，若保留小数点后，若保留小数点后10位有位有效数字则结果为效数字则结果为0.6190471917。2022-3-1611例例1.3 用迭代方法求解方程用迭代方法求解方程解：（1）对方程变形为 ，有不同的形式，比如： ； (a) ； (b) ； (c)0123xxx)(xx123xxx321xxx2111xxx2022-3-1612（2）设定初始值为）设定初始值为1，编程迭代求解，编程迭代求解lx=1;y=1;z=1;lfor k=1:25l x=x3-x2-1;l y=(

7、y2+y+1)(1/3);l z=1+1/z+1/z2;lendlx,y,z x = -Infy = 1.8393z = 1.8393 在程序中，函数在程序中，函数x,y,z分别对应方程（分别对应方程（a）（b ）（）（c），从结果可以看出方程（），从结果可以看出方程（ a）不收敛，结果趋于负无穷大，方程不收敛，结果趋于负无穷大，方程(b)(c)收敛，结果为收敛，结果为1.8393。而且，还可证明。而且，还可证明（b）比（）比（c）收敛速度快。）收敛速度快。2022-3-16131.3 方程组的求解方法方程组的求解方法l1.3.1 线性方程组的求解线性方程组的求解l对于线性方程组对于线性方程组

8、 l可以写成矩阵的形式可以写成矩阵的形式l我们在线性代数中已经学习了线性方程组的求我们在线性代数中已经学习了线性方程组的求解方法。解方法。 11 111,1 1nnmmnnma xa xba xa xb. bAX 2022-3-1614l由线性代数的知识可知，线性方程组的解可能由线性代数的知识可知，线性方程组的解可能出现三种情形：无解、有唯一解、有无穷多组出现三种情形：无解、有唯一解、有无穷多组解。这主要取决于系数矩阵解。这主要取决于系数矩阵A 的秩与增广矩阵的秩与增广矩阵（Ab）的秩是否相等、秩与变量个数是否）的秩是否相等、秩与变量个数是否相等，具体地：相等，具体地：l若若R(A)R(Ab）

9、,则无解；则无解；l若若R(A)R(Ab）n（n为变量个数），则为变量个数），则有唯一一组解；有唯一一组解；l若若R(A)R(Ab） n，则有无穷多组解。，则有无穷多组解。l求矩阵求矩阵A的秩可以很方便的用的秩可以很方便的用MATLAB的的rank(A)函数求得。函数求得。2022-3-1615l求解线性方程组的方法大致可以划分为两类：直接消去法直接消去法、迭代数值解法迭代数值解法。l直接消去法在线性代数中已经学过，这里不再赘述。l线性方程组可以看成是非线性方程组的特例，其迭代数值解法相同，在非线性方程组的迭代解法中介绍方法。2022-3-1616132 非线性方程组的迭代解法非线性方程组的迭

10、代解法l非线性方程组的一般形式为 l可以改写为等价的方程组 l用这个方程组进行迭代求得精确解。 0),(, 0),(111nnnxxfxxf),(),(1111nnnnxxgxxxgx2022-3-1617l例14 求解方程组l解：对方程组进行变形，构造如下的迭代函数：l或081008102122122121xxxxxxx8 . 01 . 01 . 08 . 01 . 01 . 01221222211xxxxxxx211222218101810 xxxxxx2022-3-1618求解方程组迭代产生的序列是数组：求解方程组迭代产生的序列是数组：l迭代程序如下：（初始值是（2,3），迭代次数为20

11、次）lx=2,3;y=2,3;lfor k=1:20la=0.1*x(1)2+0.1*x(2)2+0.8;lx(2)=0.1*x(1)*x(2)2+0.1*x(1)+0.8;lx(1)=a;lb=(10*y(2)-8)/(y(2)2+1);ly(2)=sqrt(10*y(1)-8-y(1)2);ly(1)= b;lendlx,y),( ,),(),2122211211nnxxxxxx（x = 1.0000 1.0000y = 2.1934 3.0205 2022-3-16191.3 MATLAB软件直接求解法软件直接求解法l1.3.1 任意函数方程与线性方程组可用命令任意函数方程与线性方程组可

12、用命令solve( )求解求解)( ,)( ,)(),(, 0)(, 0)(2111xfxfxfsolvexxxxfxfnnn2022-3-16201.3 MATLAB软件直接求解法软件直接求解法lsolve( )语句的用法：语句的用法：l1. 单变量方程：单变量方程：f(x)=0 l（1）符号解：）符号解：l例例13 求解方程：求解方程：l解：输入：x=solve(a*x2+b*x+c) l输出为：lx =l 1/2/a*(-b+(b2-4*a*c)(1/2)l 1/2/a*(-b-(b2-4*a*c)(1/2) 02cbxax2022-3-1621l（2）数字解：l如果不能求得精确的符号解

13、，可以计算可变精度的数值解。l例例1.4 解方程： l解：s=solve(x3-2*x2=x-1) ls = 1/6*(28+84*i*3(1/2)(1/3)+14/3/(28+84*i*3(1/2)(1/3)+2/3l -1/12*(28+84*i*3(1/2)(1/3)-7/3/(28+84*i*3(1/2)(1/3)+2/3+1/2*i*3(1/2)*(1/6*(28+84*i*3(1/2)(1/3)-14/3/(28+84*i*3(1/2)(1/3)l -1/12*(28+84*i*3(1/2)(1/3)-7/3/(28+84*i*3(1/2)(1/3)+2/3-1/2*i*3(1/2

14、)*(1/6*(28+84*i*3(1/2)(1/3)-14/3/(28+84*i*3(1/2)(1/3) 1223xxx2022-3-1622ldouble(s) lans =l 2.2470 l -0.8019 + 0.0000il 0.5550 - 0.0000i lvpa(s,10) lans =l2.246979604+.1e-9*il -.8019377358-.1866025404e-9*il .5549581322-.1339745960e-10*i l该方程无实根。2022-3-1623l（3）无穷解l例1.5 求解方程：l解：输入：solve(tan(x)-sin(x)=0

15、) l输出为：lans =l0 l该方程有无穷多个解，不能给出全部解，这里只得到其中的一个。2022-3-16242多变量方程组：多变量方程组：l例例1.6 解方程组 l解：输入：lx,y=solve(x2*y2,x-y/2-b) 0)(, 0)(1xfxfnbyxyx2022x = 0 0 b by = -2*b -2*b 0 0 v = 0, -2*b 0, -2*b b, 0 b, 0 v=x,y2022-3-16251.3.2 非线性方程组非线性方程组l非线性方程组仍然可以用solve( )求解（如上例1.6），一般给出的是数值解。也可以用fsolve( )函数求解，格式是：l这里 为

16、变量的初始值；fun为m文件的文件名。m文件如下所示：110( )0,( )0,( ,)(,)nnf xfxxxxxfsolve fun x0 x1( )(1)( );( )( )nfunctionffun xff xf nfx2022-3-1626l例例1.7 求解方程组：l解：这是一个非线性方程组。l（1）首先建立关于该方程组的M文件nxxflfunction eq=nxxf(x)lglobal number;lnumber=number+1;leq(1)=sin(x(1)+x(2)2+log(x(3)-7;leq(2)=3*x(1)+2x(2)-x(3)3+1;leq(3)=x(1)+x

17、(2)+x(3)-5;05012307lnsin32zyxzxzyxy2022-3-1627l（2）执行以下程序lglobal number;lnumber=0;lx=fsolve(nxxf,1,1,1)lnumber lx =l 0.5991 2.3959 2.0050lnumber =l 29 l这里迭代步骤为29次。2022-3-16281.3.3 多项式方程多项式方程l求解多项式方程除了可以用上面讲述的solve( ) 函数外，还可以直接用求解多项式方程的函数roots(p)。 l对于多项式方程l用： 求解l该方法可以求出方程的全部根（包含重根）。00111axaxaxammmm)(;

18、,011prootsaaaapmm2022-3-1629l例例18 求解多项式方程：l解：p=1,1,0,0,0,0,0,0,0,1;lroots(p)l运行结果：运行结果： lans =l -1.2131 l -0.9017 + 0.5753il -0.9017 - 0.5753il -0.2694 + 0.9406il -0.2694 - 0.9406il 0.4168 + 0.8419il 0.4168 - 0.8419il 0.8608 + 0.3344il 0.8608 - 0.3344i 0189 xx2022-3-16301.3.4 线性方程组线性方程组l线性方程组除了可以使用s

19、olve( ) 求解外，还可以使用其他的MATLAB命令。l将线性方程组写成矩阵形式：AX=b，就可以考虑用以下几种形式之一求解。llinsolve(A,b)；6.5版无此函数lsym(A)sym(b)；lAb；linv(A)*b； inv(A)表示A的逆矩阵，这里A必须为方阵且可逆；lpinv(A)*b； pinv(A) 表示A的逆矩阵，A可以为任意矩阵。2022-3-1631l例例1.9 求解线性方程组：AXbl(1)lA=4 1 0 ;1 -1 5;2 2 -3;b=6;14;-3;lx=Ab lx =l 1l 2l 3 2022-3-16321.4 案例详解：山崖高度案例详解：山崖高度

20、l某人站在山崖顶且身上带着一只具有跑表功能的手表，出于好奇心想用扔下一块石头听回声的方法来估计山崖的高度，假定他准确地测定时间t=4s，那么应该怎样来推算山崖的高度呢？l模型一：模型一：假定空气阻力不计，可以直接利用自由落体运动的公式来计算。l （1） l取g=9.8m/s2，则可求得 。这样计算非常简单，但略去了一些影响因素，误差较大。221gth mh4 .782022-3-1633l模型二：模型二：l除去地球吸引力外，对石块下落影响最大的当属空气阻力。根据流体力学知识，可设空气阻力正比于石块下落的速度，阻力系数K为常数，因而，由牛顿第二定律可得：l令 ，解得l代入初始条件： ，得l对v再

21、积分一次，得KvmgdtdvmFmKk/kgecvkt10)0(vkgc/122cekgtkghkt2022-3-1634l代入初始条件： ，l得l得到计算山崖高度的公式：l （2）l若设k=0.05，将t=4s代入上式，则可求得0)0(h22/kgc22kgekgtkghktmh4 .732022-3-1635l模型三：模型三：l进一步考虑，听到回声按下跑表，进一步考虑，听到回声按下跑表，t=4s包含了包含了反应时间，不妨假设平均反应时间为反应时间，不妨假设平均反应时间为0.1s，则，则扣除反应时间后扣除反应时间后下下落时间为落时间为3.9s，代入式（，代入式（2）求得求得 。mh9 .69

22、2022-3-1636l模型四：模型四：l再深入考虑回声传回来所需要的时间，为此，再深入考虑回声传回来所需要的时间，为此，令石块下落的的真正时间为令石块下落的的真正时间为 ,声音传回来的声音传回来的时间记为时间记为 ，我们知道声音在空气中的传播速，我们知道声音在空气中的传播速度是度是340m/s,因此得到一个方程组因此得到一个方程组1t2t9 . 33402122211ttthkgekgtkghkt2022-3-1637l利用MATLAB的fsolve()函数求解，先建立m文件mx4.mlfunction f=mx4(x)lglobal number;lnumber=number+1;lg=9

23、.8;lk=0.05; lf(1)=x(1)-g*x(2)/k-g*exp(-k*x(2)/k2+g/k2;lf(2)=x(1)-340*x(3);lf(3)=x(2)+x(3)-3.9;2022-3-1638l再执行以下程序：lglobal number;lnumber=0;lx=fsolve(mx4,0,0,0)lnumber l得到结果：lx = 63.5616 3.7131 0.1869lnumber = 125l即迭代125次后求出， l这样计算的过程实际上就是一个数学建模并求解的过程，通过循序渐进的计算得到了越来越精确的解。ststmh1869. 0,7131. 3,5616.63

24、212022-3-1639范例：波音公司飞机最佳定价策略范例：波音公司飞机最佳定价策略l问题问题l全球最大的飞机制造商波音公司自1955年推出的波音707开始，成功地开发了一系列的喷气式客机。问题：讨论该公司对一种新型客机最优定价策略的数学模型。l问题分析问题分析l定价策略涉及到诸多因素，这里考虑以下主要因素： l价格、竞争对手的行为、出售客机的数量、波音公司的客机制造量、制造成本、波音公司的市场占有率等等因素。2022-3-1640l假设及模型假设及模型l价格记为p，根据实际情况，对于民航飞机制造商，能够与波音公司抗衡的竞争对手只有一个，因此他们可以在价格上达成一致，具体假设如下： 1. 型

25、号：为了研究方便，假设只有一种型号飞机；2. 销售量：其销售量只受飞机价格p的影响。预测以此价格出售，该型号飞机全球销售量为N。N应该受到诸多因素的影响，假设其中价格是最主要的因素。根据市场历史的销售规律和需求曲线，假设该公司销售部门预测得到：12565578)2pppNN（2022-3-16413. 市场占有率：既然在价格上达成一致，即价格的变化是同步的，因此，不同定价不会影响波音公司的市场占有率，因此市场占有率是常数市场占有率是常数，记为h。4. 制造数量：假设制造量等于销售量，记为x。既然可以预测该型号飞机全球销售量，结合波音公司的市场占有率，可以得到)(*pNhx 2022-3-164

26、25. 制造成本：根据波音产品分析部门的估计，制造成本为：6. 利润：假设利润等于销售收入去掉成本，并且公司的最优策略原则为利润R(p)最大。4385 . 150)(xxxC)()(xCpxpR2022-3-1643l由以上简化的分析及假设得到波音公司飞机最佳定价策略的数学模型如下：l其中：l p，x，N0。 )()(xCpxpMaxR12565578)2pppNN（)(*pNhx 4385 . 150)(xxxC2022-3-1644模型求解模型求解l我们采用图形放大图形放大的方法求解。具体用Matlab 作出目标函数曲线图，得到一个直观的印象：最优定价策略下价格p大致在6到7之间；再用图形放大方法，进一步估计出 lp 6.2859，R=1780.8336 6.2856.28526.28546.28566.28586.2866.28626.28646.28666.28686.2871780.83321780.83331780.83331780.83341780.83341780.83351780.83351780.83361780.8336价 格 p利 润 曲 线 R(p)2022-3-1645lMatlab 程序如下(作函数曲线图的基本程序)lh=0.5;a=6.2