有无穷多组解。
求矩阵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.该模型本身是一个最值问题,由高等数学的知识,可以利用导数求驻点然后求最值。