一张建模方程组.docx
- 文档编号:11326044
- 上传时间:2023-02-26
- 格式:DOCX
- 页数:28
- 大小:159.99KB
一张建模方程组.docx
《一张建模方程组.docx》由会员分享,可在线阅读,更多相关《一张建模方程组.docx(28页珍藏版)》请在冰豆网上搜索。
一张建模方程组
第一章方程(组)模型
本章学习目的:
1.复习求解方程的基本原理和方法,掌握解方程的图形放大法和迭代算法;
2.能利用MATLAB软件编写迭代算法程序,了解迭代过程的图形表示;
3.熟练掌握用MATLAB软件的函数来求解方程和方程组;
4.通过范例展现求解实际问题的初步建模过程和MATLAB程序设计。
1.1引言
“方程是很多工程和科学工作的发动机”。
研究大型的土建结构、机械结构、输电网络、管道网络,研究经济规划、人口增长、种群繁殖等问题时,简单的分析可以直接归结为线性或非线性方程组,复杂一些要用到(偏)微分方程,求数值解时将转化为n非常大的方程组。
若干世纪以来,工程师和科学家花了大量的时间用于求解方程(组),数学家研究各种各样的方程求解方法。
本章我们就是要学习求解线性方程组、非线性方程(组)的方法,以及利用数学软件利用计算机对方程和方程组进行求解。
1.2方程的求解方法
考虑求方程f(x)=0的解,我们通常采用这样的几种方法:
因式分解法、图形放大法、数值迭代逼近法。
1.因式分解法:
这是我们最熟悉、常用的一种方法,这个方法的关键在分解因式,包括对多项式函数、三角函数和指数函数等的分解。
但对于无法进行分解的函数则无能为力。
2.图形放大法:
由于计算机的广泛应用,可以非常方便地作出函数f(x)的图形(曲线),找出曲线与x轴的交点的横坐标值,就可求出f(x)=0的近似根。
这些值尽管不精确,但是直观,方程有多少个根、在什么范围,一目了然。
并且可以借助于计算机使用图形局部放大功能,将根定位得更加准确。
3.数值迭代逼近法:
利用图形的方法或连续函数的零点存在性定理,可以推知f(x)在某一区间内有根,我们就可以用数值方法来求方程的根,这就是迭代逼近法。
迭代逼近法分为区间的迭代和点的迭代。
区间迭代又分为对分法和黄金分割法;点的迭代又分为简单迭代法、单点割线法、两点割线法、牛顿法等。
迭代失败后又可以采用加速迭代收敛方法。
1.2.1图形放大法
用图形放大法求解方程f(x)=0的步骤:
(1)建立坐标系,作曲线f(x);
(2)观察f(x)与x轴的交点;
(3)将其中的一个交点进行局部放大;
(4)该交点的横坐标值就是方程的一个根;
(5)对所有的交点进行相同的处理,就得到方程的所有解。
例1.1求方程
所有的根及大致分布范围,欲寻求其中的一个实根,并且达到一定的精度。
(1)画出
的图形;
x=-6:
0.01:
6;
y=x.^5+2*x.^2+4;
plot(x,y)
gridon;%画坐标格
我们可以看出方程在-2~+2范围有一个实根。
(2)逐次缩小范围得到较精确的根。
x=-2:
0.01:
2;
y=x.^5+2.*x.^2+4;
plot(x,y)
gridon
x=-2:
0.01:
-1.5;
y=x.^5+2.*x.^2+4;
plot(x,y)
gridon
x=-1.6:
0.01:
-1.5;
y=x.^5+2.*x.^2+4;
plot(x,y)
gridon
因此我们可以看出这个实根的值在-1.55~-1.5之间。
1.2.2简单迭代法
1.迭代算法步骤:
对方程f(x)=0求解
(1)对方程经过简单变形得到
(不是唯一的),x被称之为不动点;
(2)设置
为迭代初值,迭代过程为
,n=0,1,2……
(3)当两次迭代结果之差小于某个设定的误差值时,我们认为迭代结果是收敛的,可得到结果的近似值
。
例1.2求方程
的非负实根。
解:
由于函数
连续,并且在x=0和x=1处函数值符号相反,可以判断函数在区间(0,1)必有零点,即方程
在(0,1)内必然存在根。
(1)先将函数变形为
;
(2)设置迭代初值为0,编程进行迭代。
n=1;
x=0;
y=exp(x)/3;
ys=vpa(y,10);%给出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;
end
n,y,ys
n=
21
y=
0.6190
ys=
.6190471917
从该结果可以看出,迭代21次后两次迭代的结果误差值满足小于
的条件,结果收敛,迭代结果为0.6190,若保留小数点后10位有效数字则结果为0.6190471917。
例1.3用迭代方法求解方程
解:
(1)对方程变形为
,有不同的形式,比如
;(a)
;(b)
;……(c)
(2)设定初始值为1,编程迭代求解
x=1;y=1;z=1;
fork=1:
25
x=x^3-x^2-1;
y=(y^2+y+1)^(1/3);
z=1+1/z+1/z^2;
end
x,y,z
x=
-Inf
y=
1.8393
z=
1.8393
在程序中,函数x,y,z分别对应方程(a)(b)(c),从结果可以看出方程(a)不收敛,结果趋于负无穷大,方程(b)(c)收敛,结果为1.8393。
而且,还可证明(b)比(c)收敛速度快。
注:
这段程序和例1.2有所不同,这里是设定了固定的迭代次数。
2.迭代失败的改进:
加速迭代收敛方法
例1.3中方程(a)的迭代是失败的(即迭代不收敛),如何解决?
当我们遇到迭代失败时,可以采用以下的简单方法解决。
我们考虑不直接用
迭代,而用
和x的加权平均
进行迭代,
为参数。
即:
在满足
的条件下,取
,由此可以得到
。
在实际迭代过程中用
代替a,因此在迭代中
,
其加速迭代过程如下:
。
采用以上方法对例2.3中方程(a)进行改进,产生新的迭代函数
加速迭代过程为
。
再编写程序计算:
x=0;
fork=1:
20
x=(-2*x^3+x^2-1)/(-3*x^2+2*x+1);
end
x
x=
1.8393
x=100;
fork=1:
20
x=(-2*x^3+x^2-1)/(-3*x^2+2*x+1);
end
x
x=
1.8393
由此我们看到选用两个不同的初值0和100都能得到结果是1.8393,改进是非常有效的,同学们还可尝试不同的初值,观察其迭代的收敛性。
实验表明,它比(b)、(c)的收敛速度都快。
但要注意,x=1不能作为初值,这会出现被0除的错误。
通过这个例子我们看出,求解方程的迭代函数构造形式多样,不同迭代函数的收敛性和收敛速度都可能不同,需要在遇到具体问题时灵活应用。
1.3方程组的求解方法
1.3.1线性方程组的求解
我们在线性代数中已经学习了线性方程组的求解方法。
对于线性方程组
可以写成矩阵的形式
由线性代数的知识可知,线性方程组的解可能出现三种情形:
无解、有唯一解、有无穷多组解。
这主要取决于系数矩阵A的秩与增广矩阵(A︱b)的秩是否相等、秩与变量个数是否相等,具体地:
若R(A)≠R(A︱b),则
无解;
若R(A)=R(A︱b)=n(n为变量个数),则
有唯一一组解;
若R(A)=R(A︱b) 有无穷多组解。 求矩阵A的秩可以很方便的用MATLAB的rank(A)函数求得。 求解线性方程组的方法大致可以划分为两类: 直接消去法、迭代数值解法。 直接消去法在线性代数中已经学过,这里不再赘述。 线性方程组可以看成是非线性方程组的特例,其迭代数值解法相同,在非线性方程组的迭代解法中介绍方法。 1.3.2非线性方程组的迭代解法 非线性方程组的一般形式为 可以改写为等价的方程组 用这个方程组进行迭代求得精确解。 例1.4求解方程组 解: 对方程组进行变形,构造如下的迭代函数: (1) 或 (2) 思考: 迭代序列如何表示? 求解方程组迭代产生的序列是数组: 对本例选择初始点(0,0),(2,3),(8,9),…分别计算,迭代次数逐渐增加,观察结果。 迭代程序如下: (初始值是(2,3),迭代次数为20次) x=[2,3];y=[2,3]; fork=1: 20 a=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; end x,y x= 1.00001.0000 y= 2.19343.0205 从这个结果看出,x (1)=1,x (2)=1,和y (1)=2.1934,y (2)=3.0205是方程组的两组解。 问: 程序中a、b变量的作用是什么? 取消可以吗? 1.4MATLAB软件直接求解法 1.4.1任意函数方程与线性方程组可用命令solve()求解 solve()语句的用法: 1.单变量方程: f(x)=0 1)符号解: 例1.5求解方程: 解: 输入: 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))] 或输入: solve('a*x^2+b*x+c') 输出为: ans= 1/2/a*(-b+(b^2-4*a*c)^(1/2)) 1/2/a*(-b-(b^2-4*a*c)^(1/2)) 2)数字解: 如果不能求得精确的符号解,可以计算可变精度的数值解。 例1.6解方程: 解: 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)) 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 该方程有无穷多个解,不能给出全部解,这里只得到其中的一个。 2.多变量方程组: 例1.8解方程组 解: 输入: [x,y]=solve('x^2*y^2','x-y/2-b') x= 0 0 b b y= -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)首先建立关于该方程组的M文件nxxf functioneq=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; (2)执行以下程序 globalnumber; number=0; x=fsolve('nxxf',[1,1,1],1) number Optimizationterminated: first-orderoptimalityislessthanoptions.TolFun. x= 0.59912.39592.0050 number= 29 这里迭代步骤为29次。 1.4.3多项式方程 求解多项式方程除了可以用上面讲述的solve()函数外,还可以直接用求解多项式方程的函数roots(p)。 对于多项式方程 用: 求解 该方法可以求出方程的全部根(包含重根)。 例1.10求解多项式方程: 解: 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 0.8608+0.3344i 0.8608-0.3344i 当然,本例也可以用solve()函数求解。 s=solve('x^9+x^8+1') s= [.86082052816807526471321650032963+.33435225889790612035535442609416*i] [.41683400653657232118306037529804+.84191977308465856032531749536099*i] [-.26935193759657466239499735181026+.94057840102314507033098846338980*i] [-.90172773557825296622423876941201+.57531209407289260277720483797049*i] [-1.2131497230596399145540815088108] [-.90172773557825296622423876941201-.57531209407289260277720483797049*i] [-.26935193759657466239499735181026-.94057840102314507033098846338980*i] [.41683400653657232118306037529804-.84191977308465856032531749536099*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命令。 将线性方程组写成矩阵形式: AX=b,就可以考虑用以下几种形式之一求解。 (1)linsolve(A,b);6.5版无此函数 (2)sym(A)\sym(b); (3)A\b; (4)inv(A)*b;inv(A)表示A的逆矩阵,这里A必须为方阵且可逆; (5)pinv(A)*b;pinv(A)表示A的逆矩阵,A可以为任意矩阵。 例1.11求解线性方程组: AX=b (1) A=[410;1-15;22-3];b=[6;14;-3]; x=A\b x= 1 2 3 rankA=rank(A) C=[A,b]; rankAb=rank(C) rankA= 3 rankAb= 3 此例中rank(A)=rank(A|b)=3,说明方程有唯一解。 (2)A=[430;34-1;0-14];b=[24;30;-24]; x=sym(A)\sym(b) x= [3] [4] [-5] x=A\b 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) ans= 10 c=[A,b]; rank(c) ans= 10 此例中rank(A)=rank(A|b)=10,说明方程有唯一解。 1.5案例详解: 山崖高度 某人站在山崖顶且身上带着一只具有跑表功能的计算器,出于好奇心想用扔下一块石头听回声的方法来估计山崖的高度,假定他准确地测定时间t=4s,那么应该怎样来推算山崖的高度呢? 模型一: 假定空气阻力不计,可以直接利用自由落体运动的公式来计算。 (1) 取g=9.8m/s2,则可求得 。 这样计算非常简单,但略去了一些影响因素,误差较大。 模型二: 除去地球吸引力外,对石块下落影响最大的当属空气阻力。 根据流体力学知识,此时可设空气阻力正比于石块下落的速度,阻力系数K为常数,因而,由牛顿第二定律可得: 令 ,解得 代入初始条件: ,得 对v再积分一次,得 代入初始条件: ,得 得到计算山崖高度的公式: (2) 若设k=0.05,将t=4s代入式 (2),则可求得 。 模型三: 进一步考虑,听到回声按下跑表,t=4s包含了反应时间,不妨假设平均反应时间为0.1s,则扣除反应时间后落下时间为3.9s,代入式 (2)求得 。 模型四: 再深入考虑回声传回来所需要的时间,为此,令石块下落的的真正时间为 声音传回来的时间记为 ,我们知道声音在空气中的传播速度是340m/s,因此得到一个方程组 (3) 利用MATLAB的fsolve()函数求解,先建立m文件mx4.m functionf=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; 再执行以下程序: globalnumber; number=0; x=fsolve('mx4',[0,0,0]) number 得到结果: Optimizationterminated: first-orderoptimalityislessthanoptions.TolFun. x=63.56163.71310.1869 number=125 即迭代125次后求出, 。 这样计算的过程实际上就是一个数学建模并求解的过程,通过循序渐进的计算得到了越来越精确的解。 范例: 波音公司飞机最佳定价策略 全球最大的飞机制造商——波音公司自1955年推出的波音707开始,成功地开发了一系列的喷气式客机。 问题: 讨论该公司对一种新型客机最优定价策略的数学模型。 问题分析 定价策略涉及到诸多因素,这里考虑以下主要因素: 价格、竞争对手的行为、出售客机的数量、波音公司的客机制造量、制造成本、波音公司的市场占有率等等因素。 假设及模型 价格记为p,根据实际情况,对于民航飞机制造商,能够与波音公司抗衡的竞争对手只有一个,因此他们可以在价格上达成一致,具体假设如下: 1.型号: 为了研究方便,假设只有一种型号飞机; 2.销售量: 其销售量只受飞机价格p的影响。 预测以此价格出售,该型号飞机全球销售量为N。 N应该受到诸多因素的影响,假设其中价格是最主要的因素。 根据市场历史的销售规律和需求曲线,假设该公司销售部门预测得到: 3.市场占有率: 既然在价格上达成一致,即价格的变化是同步的,因此,不同定价不会影响波音公司的市场占有率,因此市场占有率是常数,记为h。 4.制造数量: 假设制造量等于销售量,记为x。 既然可以预测该型号飞机全球销售量,结合波音公司的市场占有率,可以得到 5.制造成本: 根据波音产品分析部门的估计,制造成本为: 6.利润: 假设利润等于销售收入去掉成本,并且公司的最优策略原则为利润R(p)最大。 由以上简化的分析及假设得到波音公司飞机最佳定价策略的数学模型如下: 其中: ; ; ,p,x,N≥0。 模型求解 我们采用图形放大的方法求解。 具体用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; fori=1: n+1 pr(i)=a+(i-1)*d;p=pr(i); n=-78*p^2+655*p+125;x=h*n;r=p*x; c=50+1.5*x+8*x^(3/4);l(i)=r-c; end plot(pr,l);gridon xlabel('价格p');title('利润曲线R(p)') 注意: 1.根据图形的具体情况,不断修改上面程序中的最长一条语句,就可以不断地放大图形,将最优解的范围限制得越来越小,直至找出满意的近似解。 2.以上的市场占有率h=0.5;对于市场占有率h的其它取值,可以类似地进行。 进一步思考的几个问题 1.求出h取其它值时的最优价格,并进行比较; 2.该模型本身是一个最值问题,由高等数学的知识,可以利用导数求驻点然后求最值。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 一张 建模 方程组