一张建模方程组.docx

1、一张建模方程组第一章 方程（组）模型本章学习目的：1复习求解方程的基本原理和方法，掌握解方程的图形放大法和迭代算法；2能利用MATLAB软件编写迭代算法程序，了解迭代过程的图形表示；3熟练掌握用MATLAB软件的函数来求解方程和方程组；4通过范例展现求解实际问题的初步建模过程和MATLAB程序设计。1.1 引言 “方程是很多工程和科学工作的发动机”。研究大型的土建结构、机械结构、输电网络、管道网络，研究经济规划、人口增长、种群繁殖等问题时，简单的分析可以直接归结为线性或非线性方程组，复杂一些要用到（偏）微分方程，求数值解时将转化为n非常大的方程组。若干世纪以来，工程师和科学家花了大量的时间用于

2、求解方程（组），数学家研究各种各样的方程求解方法。本章我们就是要学习求解线性方程组、非线性方程（组）的方法，以及利用数学软件利用计算机对方程和方程组进行求解。1.2 方程的求解方法 考虑求方程f(x)=0的解，我们通常采用这样的几种方法：因式分解法、图形放大法、数值迭代逼近法。1因式分解法：这是我们最熟悉、常用的一种方法，这个方法的关键在分解因式，包括对多项式函数、三角函数和指数函数等的分解。但对于无法进行分解的函数则无能为力。2图形放大法：由于计算机的广泛应用，可以非常方便地作出函数f(x)的图形（曲线），找出曲线与x轴的交点的横坐标值，就可求出f(x)0的近似根。这些值尽管不精确，但是直观

3、，方程有多少个根、在什么范围，一目了然。并且可以借助于计算机使用图形局部放大功能，将根定位得更加准确。3数值迭代逼近法：利用图形的方法或连续函数的零点存在性定理，可以推知f(x)在某一区间内有根，我们就可以用数值方法来求方程的根，这就是迭代逼近法。 迭代逼近法分为区间的迭代和点的迭代。 区间迭代又分为对分法和黄金分割法；点的迭代又分为简单迭代法、单点割线法、两点割线法、牛顿法等。迭代失败后又可以采用加速迭代收敛方法。1.2.1图形放大法用图形放大法求解方程f(x)0的步骤：（1） 建立坐标系，作曲线f(x)；（2） 观察f(x)与x轴的交点；（3） 将其中的一个交点进行局部放大；（4） 该交点

4、的横坐标值就是方程的一个根；（5） 对所有的交点进行相同的处理，就得到方程的所有解。例1.1 求方程所有的根及大致分布范围，欲寻求其中的一个实根，并且达到一定的精度。（1） 画出的图形；x=-6:0.01:6;y=x.5+2*x.2+4;plot(x,y) grid on; %画坐标格 我们可以看出方程在22范围有一个实根。（2） 逐次缩小范围得到较精确的根。x=-2:0.01:2;y=x.5+2.*x.2+4;plot(x,y)grid on x=-2:0.01:-1.5;y=x.5+2.*x.2+4;plot(x,y)grid on x=-1.6:0.01:-1.5;y=x.5+2.*x.

5、2+4;plot(x,y) grid on 因此我们可以看出这个实根的值在-1.55-1.5之间。1.2.2 简单迭代法1迭代算法步骤： 对方程f(x)=0求解（1） 对方程经过简单变形得到（不是唯一的），x 被称之为不动点；（2） 设置为迭代初值，迭代过程为，n=0,1,2（3） 当两次迭代结果之差小于某个设定的误差值时，我们认为迭代结果是收敛的，可得到结果的近似值。例12 求方程的非负实根。解：由于函数连续，并且在x=0和x=1处函数值符号相反，可以判断函数在区间（0，1）必有零点，即方程在（0，1）内必然存在根。（1） 先将函数变形为；（2） 设置迭代初值为0，编程进行迭代。n=1;x=

6、0;y=exp(x)/3;ys=vpa(y,10); %给出y的数值型结果，有效位数为10z=abs(y-x);while z10(-5) x=y; y=exp(x)/3; ys=vpa(y,10); z=abs(y-x); n=n+1;endn,y,ys n = 21y = 0.6190ys =.6190471917 从该结果可以看出，迭代21次后两次迭代的结果误差值满足小于的条件，结果收敛，迭代结果为0.6190，若保留小数点后10位有效数字则结果为0.6190471917。例1.3 用迭代方法求解方程解：（1）对方程变形为，有不同的形式，比如； (a)； (b)； (c)（2）设定初始值

7、为1，编程迭代求解x=1;y=1;z=1;for k=1:25x=x3-x2-1;y=(y2+y+1)(1/3);z=1+1/z+1/z2;endx,y,z x = -Infy = 1.8393z = 1.8393 在程序中，函数x,y,z分别对应方程（a）（b ）（c），从结果可以看出方程（ a）不收敛，结果趋于负无穷大，方程(b)(c)收敛，结果为1.8393。而且，还可证明（b）比（c）收敛速度快。注：这段程序和例1.2有所不同，这里是设定了固定的迭代次数。2迭代失败的改进：加速迭代收敛方法例1.3中方程（a）的迭代是失败的（即迭代不收敛），如何解决？当我们遇到迭代失败时，可以采用以下的

8、简单方法解决。我们考虑不直接用迭代，而用和x的加权平均进行迭代，为参数。即： 在满足的条件下，取，由此可以得到。在实际迭代过程中用代替a，因此在迭代中，其加速迭代过程如下：。采用以上方法对例2.3中方程（a）进行改进，产生新的迭代函数 加速迭代过程为。再编写程序计算：x=0;for k=1:20x=(-2*x3+x2-1)/(-3*x2+2*x+1);endx x = 1.8393 x=100;for k=1:20x=(-2*x3+x2-1)/(-3*x2+2*x+1);endx x = 1.8393 由此我们看到选用两个不同的初值0和100都能得到结果是1.8393，改进是非常有效的，同学们

9、还可尝试不同的初值，观察其迭代的收敛性。实验表明，它比（b）、（c）的收敛速度都快。但要注意，x=1不能作为初值，这会出现被0除的错误。通过这个例子我们看出，求解方程的迭代函数构造形式多样，不同迭代函数的收敛性和收敛速度都可能不同，需要在遇到具体问题时灵活应用。1.3 方程组的求解方法1.3.1 线性方程组的求解 我们在线性代数中已经学习了线性方程组的求解方法。对于线性方程组可以写成矩阵的形式由线性代数的知识可知，线性方程组的解可能出现三种情形：无解、有唯一解、有无穷多组解。这主要取决于系数矩阵A 的秩与增广矩阵（Ab）的秩是否相等、秩与变量个数是否相等，具体地：若R(A)R(Ab）,则无解；

10、若R(A)R(Ab）n（n为变量个数），则有唯一一组解；若R(A)R(Ab） n，则有无穷多组解。求矩阵A的秩可以很方便的用MATLAB的rank(A)函数求得。求解线性方程组的方法大致可以划分为两类：直接消去法、迭代数值解法。直接消去法在线性代数中已经学过，这里不再赘述。线性方程组可以看成是非线性方程组的特例，其迭代数值解法相同，在非线性方程组的迭代解法中介绍方法。132 非线性方程组的迭代解法非线性方程组的一般形式为可以改写为等价的方程组用这个方程组进行迭代求得精确解。例14 求解方程组解：对方程组进行变形，构造如下的迭代函数： （1）或 （2）思考：迭代序列如何表示？求解方程组迭代产生的

11、序列是数组： 对本例选择初始点（0，0），（2，3），（8，9），分别计算，迭代次数逐渐增加，观察结果。迭代程序如下：（初始值是（2,3），迭代次数为20次）x=2,3;y=2,3;for k=1:20a=0.1*x(1)2+0.1*x(2)2+0.8;x(2)=0.1*x(1)*x(2)2+0.1*x(1)+0.8;x(1)=a;b=(10*y(2)-8)/(y(2)2+1);y(2)=sqrt(10*y(1)-8-y(1)2);y(1)= b;endx,y x = 1.0000 1.0000y = 2.1934 3.0205 从这个结果看出，x(1)=1,x(2)=1,和y(1)=2.19

12、34,y(2)=3.0205是方程组的两组解。 问：程序中a、b变量的作用是什么？取消可以吗？1.4 MATLAB软件直接求解法1.4.1 任意函数方程与线性方程组可用命令solve( )求解solve( )语句的用法：1 单变量方程：f(x)=01） 符号解：例15 求解方程： 解：输入：x=solve(a*x2+b*x+c) 输出为：x = 1/2/a*(-b+(b2-4*a*c)(1/2) 1/2/a*(-b-(b2-4*a*c)(1/2) 或输入：solve(a*x2+b*x+c) 输出为：ans = 1/2/a*(-b+(b2-4*a*c)(1/2) 1/2/a*(-b-(b2-4*

13、a*c)(1/2)2） 数字解：如果不能求得精确的符号解，可以计算可变精度的数值解。例1.6 解方程： 解：s=solve(x3-2*x2=x-1) s = 1/6*(28+84*i*3(1/2)(1/3)+14/3/(28+84*i*3(1/2)(1/3)+2/3 -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) -1/12*(28+84*i*3(1/2)(1/3)-7/3/(28+84*i*3(1/2

14、)(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) double(s) ans = 2.2470 -0.8019 + 0.0000i 0.5550 - 0.0000i vpa(s,10) ans = 2.246979604+.1e-9*i -.8019377358-.1866025404e-9*i .5549581322-.1339745960e-10*i 该方程无实根。 3）无穷解 例1.7 求解方程： 解：输入：solve(tan(x)-sin(x)=0) 输出为：ans =0 该方程有无穷

15、多个解，不能给出全部解，这里只得到其中的一个。2多变量方程组： 例1.8 解方程组解：输入：x,y=solve(x2*y2,x-y/2-b) x = 0 0 b by = -2*b -2*b 0 0 v=x,y v = 0, -2*b 0, -2*b b, 0 b, 0 1.4.2 非线性方程组 非线性方程组仍然可以用solve( )求解（如上例1.8），一般给出的是数值解。也可以用fsolve( )函数求解，格式是：这里为变量的初始值；options可缺省，若options=1，表示输出中间结果；fun为m文件的文件名。m文件如下所示：例1.9 求解方程组： 解：这是一个非线性方程组。（1）

16、首先建立关于该方程组的M文件nxxffunction eq=nxxf(x)global number;number=number+1;eq(1)=sin(x(1)+x(2)2+log(x(3)-7;eq(2)=3*x(1)+2x(2)-x(3)3+1;eq(3)=x(1)+x(2)+x(3)-5;（2）执行以下程序global number;number=0;x=fsolve(nxxf,1,1,1,1)number Optimization terminated: first-order optimality is less than options.TolFun.x = 0.5991 2.3

17、959 2.0050number = 29 这里迭代步骤为29次。1.4.3 多项式方程 求解多项式方程除了可以用上面讲述的solve( ) 函数外，还可以直接用求解多项式方程的函数roots(p)。 对于多项式方程用： 求解该方法可以求出方程的全部根（包含重根）。例110 求解多项式方程： 解：p=1,1,0,0,0,0,0,0,0,1 ;roots(p) ans = -1.2131 -0.9017 + 0.5753i -0.9017 - 0.5753i -0.2694 + 0.9406i -0.2694 - 0.9406i 0.4168 + 0.8419i 0.4168 - 0.8419i

18、 0.8608 + 0.3344i 0.8608 - 0.3344i 当然，本例也可以用solve( ) 函数求解。s=solve(x9+x8+1) s = .86082052816807526471321650032963+.33435225889790612035535442609416*i .41683400653657232118306037529804+.84191977308465856032531749536099*i -.26935193759657466239499735181026+.94057840102314507033098846338980*i -.90172773

19、557825296622423876941201+.57531209407289260277720483797049*i -1.2131497230596399145540815088108 -.90172773557825296622423876941201-.57531209407289260277720483797049*i -.26935193759657466239499735181026-.94057840102314507033098846338980*i .41683400653657232118306037529804-.841919773084658560325317495

20、36099*i .86082052816807526471321650032963-.33435225889790612035535442609416*i double(s) ans = 0.8608 + 0.3344i 0.4168 + 0.8419i -0.2694 + 0.9406i -0.9017 + 0.5753i -1.2131 -0.9017 - 0.5753i -0.2694 - 0.9406i 0.4168 - 0.8419i 0.8608 - 0.3344i 1.4.4 线性方程组 线性方程组除了可以使用solve( ) 求解外，还可以使用其他的MATLAB命令。 将线性方

21、程组写成矩阵形式：AX=b，就可以考虑用以下几种形式之一求解。（1） linsolve(A,b)；6.5版无此函数（2） sym(A)sym(b)；（3） Ab；（4） inv(A)*b； inv(A)表示A的逆矩阵，这里A必须为方阵且可逆；（5） pinv(A)*b； pinv(A) 表示A的逆矩阵，A可以为任意矩阵。例1.11 求解线性方程组：AXb(1)A=4 1 0 ;1 -1 5;2 2 -3;b=6;14;-3;x=Ab x = 1 2 3 rankA=rank(A)C=A,b;rankAb=rank(C) rankA = 3rankAb = 3 此例中rank(A)=rank(A

22、|b)=3，说明方程有唯一解。 (2) A=4 3 0;3 4 -1;0 -1 4;b=24;30;-24;x=sym(A)sym(b) x = 3 4 -5 x=Ab x = 3.0000 4.0000 -5.0000 rank(A) ans = 3 c=A,b; rank(c) ans = 3 此例中rank(A)=rank(A|b)=3，说明方程有唯一解。(3) A=hilb(10);b=sum(A); x=inv(A)*b x = 1.0000 1.0000 1.0000 1.0000 0.9999 1.0003 0.9996 1.0009 0.9998 1.0000 rank(A)

23、ans = 10 c=A,b;rank(c) ans = 10 此例中rank(A)=rank(A|b)=10，说明方程有唯一解。1.5 案例详解：山崖高度某人站在山崖顶且身上带着一只具有跑表功能的计算器，出于好奇心想用扔下一块石头听回声的方法来估计山崖的高度，假定他准确地测定时间t=4s，那么应该怎样来推算山崖的高度呢？模型一：假定空气阻力不计，可以直接利用自由落体运动的公式来计算。 （1）取g=9.8m/s2，则可求得。这样计算非常简单，但略去了一些影响因素，误差较大。模型二：除去地球吸引力外，对石块下落影响最大的当属空气阻力。根据流体力学知识，此时可设空气阻力正比于石块下落的速度，阻力系

24、数K为常数，因而，由牛顿第二定律可得：令，解得代入初始条件：，得对v再积分一次，得代入初始条件：，得得到计算山崖高度的公式： （2）若设k=0.05，将t=4s代入式（2），则可求得。模型三：进一步考虑，听到回声按下跑表，t=4s包含了反应时间，不妨假设平均反应时间为0.1s，则扣除反应时间后落下时间为3.9s，代入式（2）求得。模型四：再深入考虑回声传回来所需要的时间，为此，令石块下落的的真正时间为,声音传回来的时间记为，我们知道声音在空气中的传播速度是340m/s,因此得到一个方程组 （3）利用MATLAB的fsolve()函数求解，先建立m文件mx4.mfunction f=mx4(x)

25、global number;number=number+1;g=9.8;k=0.05; f(1)=x(1)-g*x(2)/k-g*exp(-k*x(2)/k2+g/k2;f(2)=x(1)-340*x(3);f(3)=x(2)+x(3)-3.9;再执行以下程序：global number;number=0;x=fsolve(mx4,0,0,0)number 得到结果：Optimization terminated: first-order optimality is less than options.TolFun.x = 63.5616 3.7131 0.1869number = 125即迭

26、代125次后求出，。这样计算的过程实际上就是一个数学建模并求解的过程，通过循序渐进的计算得到了越来越精确的解。范例：波音公司飞机最佳定价策略 全球最大的飞机制造商波音公司自1955年推出的波音707开始，成功地开发了一系列的喷气式客机。问题：讨论该公司对一种新型客机最优定价策略的数学模型。 问题分析定价策略涉及到诸多因素，这里考虑以下主要因素： 价格、竞争对手的行为、出售客机的数量、波音公司的客机制造量、制造成本、波音公司的市场占有率等等因素。假设及模型价格记为p，根据实际情况，对于民航飞机制造商，能够与波音公司抗衡的竞争对手只有一个，因此他们可以在价格上达成一致，具体假设如下： 1. 型号：

27、为了研究方便，假设只有一种型号飞机；2. 销售量：其销售量只受飞机价格p的影响。预测以此价格出售，该型号飞机全球销售量为N。N应该受到诸多因素的影响，假设其中价格是最主要的因素。根据市场历史的销售规律和需求曲线，假设该公司销售部门预测得到：3. 市场占有率：既然在价格上达成一致，即价格的变化是同步的，因此，不同定价不会影响波音公司的市场占有率，因此市场占有率是常数，记为h。4. 制造数量：假设制造量等于销售量，记为x。既然可以预测该型号飞机全球销售量，结合波音公司的市场占有率，可以得到5. 制造成本：根据波音产品分析部门的估计，制造成本为：6. 利润：假设利润等于销售收入去掉成本，并且公司的最

28、优策略原则为利润R(p)最大。由以上简化的分析及假设得到波音公司飞机最佳定价策略的数学模型如下：其中：；，p，x，N0。模型求解我们采用图形放大的方法求解。具体用Matlab 作出目标函数曲线图，得到一个直观的印象：最优定价策略下价格p大致在6到7之间；再用图形放大方法，进一步估计出(如下图)p 6.2859，R=1780.8336。Matlab 程序如下(作函数曲线图的基本程序)h=0.5;a=6.285;b=6.287;n=80;d=(b-a)/n;for i=1:n+1pr(i)=a+(i-1)*d;p=pr(i);n=-78*p2+655*p+125;x=h*n;r=p*x;c=50+1.5*x+8*x(3/4);l(i)=r-c;endplot(pr,l);grid onxlabel(价格p);title(利润曲线R(p) 注意：1. 根据图形的具体情况，不断修改上面程序中的最长一条语句，就可以不断地放大图形，将最优解的范围限制得越来越小，直至找出满意的近似解。 2. 以上的市场占有率h=0.5；对于市场占有率h的其它取值，可以类似地进行。进一步思考的几个问题 1 求出h取其它值时的最优价格，并进行比较； 2 该模型本身是一个最值问题，由高等数学的知识，可以利用导数求驻点然后求最值。