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

第一章方程（组）模型2023/2/41本章学习目的：本章我们就是要学习求解线性方程组、非线性方程（组）的方法，以及利用数学软件和计算机对方程和方程组进行求解。考虑求方程f(x)=0的解，我们通常采用这样的几种方法：因式分解法、图形放大法、数值迭代逼近法。2023/2/42本章学习目的：复习求解方程的基本原理和方法，掌握解方程的图形放大法和迭代算法；熟练掌握用MATLAB软件的函数来求解方程和方程组；通过范例展现求解实际问题的初步建模过程和MATLAB程序设计。2023/2/431.1图形放大法：由于计算机的广泛应用，可以非常方便地作出函数f(x)的图形（曲线），找出曲线与x轴的交点的横坐标值，就可求出f(x)＝0的近似根。这些值尽管不精确，但是直观，方程有多少个根、在什么范围，一目了然。并且可以借助于计算机使用图形局部放大功能，将根定位得更加准确。2023/2/44用图形放大法求解方程f(x)＝0的步骤：建立坐标系，作曲线f(x)；观察f(x)与x轴的交点；将其中的一个交点进行局部放大；该交点的横坐标值就是方程的一个根；对所有的交点进行相同的处理，就得到方程的所有解。2023/2/45例1.1求方程所有的根及大致分布范围，欲寻求其中的一个实根，并且达到一定的精度。画出的图形；程序如下：x=-6:0.01:6;y=x.^5+2*x.^2+4;plot(x,y)gridon;%画坐标格2023/2/46我们可以看出方程在－2～＋2范围有一个实根。2023/2/47逐次缩小范围得到较精确的根。最后我们可以看出这个实根的值在-1.56~-1.54之间。2023/2/481.2简单迭代法迭代算法步骤：对方程f(x)=0求解对方程经过简单变形得到（不是唯一的），x被称之为不动点；设置为迭代初值，迭代过程为，n=0,1,2……当两次迭代结果之差小于某个设定的误差值时，我们认为迭代结果是收敛的，可得到结果的近似值。2023/2/49例1.2求方程的非负实根。解：由于函数连续，并且在x=0和x=1处函数值符号相反，可以判断函数在区间（0，1）必有零点，即方程在（0，1）内必然存在根。先将函数变形为；设置迭代初值为0，编程进行迭代。2023/2/410n=1;x=0;y=exp(x)/3;ys=vpa(y,10);%ys取y的10位有效数字z=abs(y-x);whilez>10^(-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。2023/2/411例1.3用迭代方法求解方程解：（1）对方程变形为，有不同的形式，比如：；(a)

；(b)

；……(c)2023/2/412（2）设定初始值为1，编程迭代求解x=1;y=1;z=1;fork=1:25x=x^3-x^2-1;y=(y^2+y+1)^(1/3);z=1+1/z+1/z^2;endx,y,z

x=-Infy=1.8393z=1.8393在程序中，函数x,y,z分别对应方程（a）（b）（c），从结果可以看出方程（a）不收敛，结果趋于负无穷大，方程(b)(c)收敛，结果为1.8393。而且，还可证明（b）比（c）收敛速度快。2023/2/4131.3方程组的求解方法1.3.1线性方程组的求解对于线性方程组

可以写成矩阵的形式我们在线性代数中已经学习了线性方程组的求解方法。

2023/2/414由线性代数的知识可知，线性方程组的解可能出现三种情形：无解、有唯一解、有无穷多组解。这主要取决于系数矩阵A的秩与增广矩阵（A︱b）的秩是否相等、秩与变量个数是否相等，具体地：若R(A)≠R(A︱b）,则无解；若R(A)＝R(A︱b）＝n（n为变量个数），则有唯一一组解；若R(A)＝R(A︱b）<n，则有无穷多组解。求矩阵A的秩可以很方便的用MATLAB的rank(A)函数求得。2023/2/415求解线性方程组的方法大致可以划分为两类：直接消去法、迭代数值解法。直接消去法在线性代数中已经学过，这里不再赘述。线性方程组可以看成是非线性方程组的特例，其迭代数值解法相同，在非线性方程组的迭代解法中介绍方法。2023/2/4161．3．2非线性方程组的迭代解法非线性方程组的一般形式为可以改写为等价的方程组用这个方程组进行迭代求得精确解。2023/2/417例1．4求解方程组解：对方程组进行变形，构造如下的迭代函数：或2023/2/418求解方程组迭代产生的序列是数组：迭代程序如下：（初始值是（2,3），迭代次数为20次）x=[2,3];y=[2,3];fork=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.00001.0000y=2.19343.02052023/2/4191.3MATLAB软件直接求解法1.3.1任意函数方程与线性方程组可用命令solve()求解2023/2/4201.3MATLAB软件直接求解法solve()语句的用法：1.单变量方程：f(x)=0（1）符号解：例1．3求解方程：解：输入：x=solve('a*x^2+b*x+c')

输出为：x=[1/2/a*(-b+(b^2-4*a*c)^(1/2))][1/2/a*(-b-(b^2-4*a*c)^(1/2))]2023/2/421（2）数字解：如果不能求得精确的符号解，可以计算可变精度的数值解。例1.4

解方程：解：s=solve('x^3-2*x^2=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))^(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))2023/2/422double(s)

ans=2.2470-0.8019+0.0000i0.5550-0.0000ivpa(s,10)

ans=2.246979604+.1e-9*i-.8019377358-.1866025404e-9*i.5549581322-.1339745960e-10*i该方程无实根。2023/2/423（3）无穷解例1.5求解方程：解：输入：solve('tan(x)-sin(x)=0')输出为：ans=0 该方程有无穷多个解，不能给出全部解，这里只得到其中的一个。2023/2/4242．多变量方程组：例1.6

解方程组解：输入：[x,y]=solve('x^2*y^2','x-y/2-b')

x=00bby=-2*b-2*b00v=[0,-2*b][0,-2*b][b,0][b,0]v=[x,y]2023/2/4251.3.2非线性方程组非线性方程组仍然可以用solve()求解（如上例1.6），一般给出的是数值解。也可以用fsolve()函数求解，格式是：这里

为变量的初始值；fun为m文件的文件名。m文件如下所示：2023/2/426例1.7

求解方程组：解：这是一个非线性方程组。（1）首先建立关于该方程组的M文件nxxffunctioneq=nxxf(x)globalnumber;number=number+1;eq(1)=sin(x(1))+x(2)^2+log(x(3))-7;eq(2)=3*x(1)+2^x(2)-x(3)^3+1;eq(3)=x(1)+x(2)+x(3)-5;2023/2/427（2）执行以下程序globalnumber;number=0;x=fsolve('nxxf',[1,1,1])number

x=0.59912.39592.0050number=29这里迭代步骤为29次。2023/2/4281.3.3多项式方程求解多项式方程除了可以用上面讲述的solve()函数外，还可以直接用求解多项式方程的函数roots(p)。

对于多项式方程用：求解该方法可以求出方程的全部根（包含重根）。2023/2/429例1．8

求解多项式方程：解：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.9406i0.4168+0.8419i0.4168-0.8419i0.8608+0.3344i0.8608-0.3344i2023/2/4301.3.4线性方程组线性方程组除了可以使用solve()求解外，还可以使用其他的MATLAB命令。 将线性方程组写成矩阵形式：AX=b，就可以考虑用以下几种形式之一求解。linsolve(A,b)；6.5版无此函数sym(A)\sym(b)；A\b；inv(A)*b；inv(A)表示A的逆矩阵，这里A必须为方阵且可逆；pinv(A)*b；pinv(A)表示A的逆矩阵，A可以为任意矩阵。2023/2/431例1.9

求解线性方程组：AX＝b(1)A=[410;1-15;22-3];b=[6;14;-3];x=A\b

x=1232023/2/4321.4案例详解：山崖高度某人站在山崖顶且身上带着一只具有跑表功能的手表，出于好奇心想用扔下一块石头听回声的方法来估计山崖的高度，假定他准确地测定时间t=4s，那么应该怎样来推算山崖的高度呢？模型一：假定空气阻力不计，可以直接利用自由落体运动的公式来计算。（1）取g=9.8m/s2，则可求得。这样计算非常简单，但略去了一些影响因素，误差较大。2023/2/433模型二：除去地球吸引力外，对石块下落影响最大的当属空气阻力。根据流体力学知识，可设空气阻力正比于石块下落的速度，阻力系数K为常数，因而，由牛顿第二定律可得：令，解得代入初始条件：，得对v再积分一次，得2023/2/434代入初始条件：，得得到计算山崖高度的公式：（2）若设k=0.05，将t=4s代入上式，则可求得2023/2/435模型三：进一步考虑，听到回声按下跑表，t=4s包含了反应时间，不妨假设平均反应时间为0.1s，则扣除反应时间后下落时间为3.9s，代入式（2）求得。2023/2/436模型四：再深入考虑回声传回来所需要的时间，为此，令石块下落的的真正时间为,声音传回来的时间记为，我们知道声音在空气中的传播速度是340m/s,因此得到一个方程组2023/2/437利用MATLAB的fsolve()函数求解，先建立m文件mx4.mfunctionf=mx4(x)globalnumber;number=number+1;g=9.8;k=0.05;

f(1)=x(1)-g*x(2)/k-g*exp(-k*x(2))/k^2+g/k^2;f(2)=x(1)-340*x(3);f(3)=x(2)+x(3)-3.9;2023/2/438再执行以下程序：globalnumber;number=0;x=fsolve('mx4',[0,0,0])number

得到结果：x=63.56163.71310.1869number=125即迭代125次后求出，这样计算的过程实际上就是一个数学建模并求解的过程，通过循序渐进的计算得到了越来越精确的解。2023/2/439范例：波音公司飞机最佳定价策略问题全球最大的飞机制造商——波音公司自1955年推出的波音707开始，成功地开发了一系列的喷气式客机。问题：讨论该公司对一种新型客机最优定价策略的数学模型。问题分析定价策略涉及到诸多因素，这里考虑以下主要因素：价格、竞争对手的行为、出售客机的数量、波音公司的客机制造量、制造成本、波音公司的市场占有率等等因素。2023/2/440假设及模型价格记为p，根据实际情况，对于民航飞机制造商，能够与波音公司抗衡的竞争对手只有一个，因此他们可以在价格上达成一致，具体假设如下：1.型号：为了研究方便，假设只有一种型号飞机；2.销售量：其销售量只受飞机价格p的影响。预测以此价格出售，该型号飞机全球销售量为N。N应该受到诸多因素的影响，假设其中价格是最主要的因素。根据市场历史的销售规律和需求曲线，假设该公司销售部门预测得到：2023/2/4413.市场占有率：既然在价格上达成一致，即价格的变化是同步的，因此，不同定价不会影响波音公司的市场占有率，因此市场占有率是常数，记为h。4.制造数量：假设制造量等于销售量，记为x。既然可以预测该型号飞机全球销售量，结合波音公司的市场占有率，可以得到2023/2/4425.制造成本：根据波音产品分析部门的估计，制造成本为：6.利润：假设利润等于销售收入去掉成本，并且公司的最优策略原则为利润R(p)最大。2023/2/443由以上简化的分析及假设得到波音公司飞机最佳定价策略的数学模型如下：其中：

p，x，N≥0。2023/2/444模型求解我们采用图形放大的方法求解。具体用Matlab作出目标函数曲线图，得到一个直观的印象：最优定价策略下价格p大致在6到7之间；再用图形放大方法，进一步估计出p≈6.2859，R=1780.83362023/2/445Matlab程序如下(作函数曲线图的基本程序)h=0.5;a=6.285;b=6.287;n=80;d=(b-a)/n;fori=1:n+1pr(i)=a+(i-1)*d;p=pr(i);n=-78*p^2+655*p+125;x=h*n;r=p*x;c=50+