七 代数方程数值解

1x?4x?y?e?1??10例1：求方程组?的实根，

??x?4y?1x2?0?8?syms x y;

s=solve(4*x-y+exp(x)/10-1,-x+4*y+x^2/8,x,y);

结果为复数，但可验证x=0.2326,y=0.0565,是方程组的解。

求方程的数值解用命令：fzero, fsolve

fzero的应用格式为x=fzero(@fun,x0)，fun是一元函数，输出x是函数fun的一个零点。x0是一个数时，输出x0附近的零点；x0是向量[a,b]时，输出函数在区间[a,b]中的零点。

例2：求函数y?xsin(x2?x?1)在-2与-0.1之间的零点。 先定义函数

function y=my(x) %函数名起为my y=x*sin(x^2-x-1);

再解方程 x=fzero(@my,0)得x=0;

再解 x=fzero(@my,-2)得x=-2.2447； 都不满足，画图，观察，解

x=fzero(@my,-0.6)得x=-0.6180; x=fzero(@my,-1.6)得x=-1.5956;

如果用x=fzero(@my,[-2,-0.1]),则出现 ??? Error using ==> fzero

The function values at the interval endpoints must differ in sign.

所以仍要通过画图观察求解。

fsolve可以解方程组，应用格式为[x,f,h]=fsolve(@fun,x0)，输出x是函数fun在x0附近的一个零点，x0为迭代初值；输出f为fun在x处的函数值，应接近于0；h大于0说明计算结果可靠，否则不可靠，需另寻他法，函数fzero也有类似输出。

例3：解例2，[x,f,h]=fsolve(@fun,-1.6)， 例4：解例1，先定义函数 function y=you(x)

y(1)=4*x(1)-x(2)+exp(x(1))/10-1; y(2)=-x(1)+4*x(2)+x(1)^2/8; 取迭代初值x0=[0,0]，

解方程组[x,f,h]=fsolve(@you,x0)， 或[x,f,h]=fsolve(@you,[0,0])

练习1：fzero和fsolve只能求实根，试用它们解x?x?1?0，看看发生了什么，再上网查函数roots的用法，用它来解方程x?x?1?0，再用函数solve求解，比较有什么不同。

22练习2：两个椭圆可能有0到4个交点，求下列两个椭圆的所有交点的坐标：

(x?2)2?(y?3?2x)2?5；2(x?3)2?(y/3)2?4

先用ezplot画图观察交点的大概位置，再用fsolve解出交点

八 数值积分

符号积分：int 例1：求

?1?11?x2dx， ?111?1?x2?1dx

方法一：I=int('1/(1+sqrt(1-x^2))')

pretty(I) 方法二：syms x;

f=1/(1+sqrt(1-x^2));

I=int(f)

数值积分：

梯形法：I=trapz(x,y)，x是积分区间向量，y是被积函数值； 抛物线法（Simpson法）： I=quad(@fun,a,b)，fun是被积函数，用法与fzero相同，

a,b是积分上下限；

高精度Lobatto法：I=quadl(@fun,a,b)，用法同quad 例2：正态分布

?1?1e1?x22dx，信号处理?sinxdx，

?1x1椭圆周长L?2以

??0a2sin2t?b2cos2tdt

?1?1e1?x22dx为例：

x=-1:0.5:1; %步长为0.5，从-0.1到0.1积分 y=exp(-x.^2/2); I=trapz(x,y)

x1=-1:0.05:1; %画图比较 y1=exp(-x1.^2/2); plot(x,y,x1,y1)

I1=trapz(x1,y1) %步长改为0.05积分 I2=int('exp(-x^2/2)',-1,1) %符号积分

vpa(I2) %符号积分的结果显示为数值，注意结果仍为符号变量

练习：我国第一颗人造地球卫星近地点距离h=439km，远地点距离H=2384km，绕地球一周需 114分钟。取地球长半轴R=6400km。由于地球位于卫星椭圆轨道的一个焦点上，根据近地点距离和远地点距离可分别计算出椭圆长半轴、椭圆半焦距、椭圆短半轴

a?(h?H?2R)/2，c?(H?h)/2，b?a2?c2，

利用椭圆周长公式：L?2??0a2sin2t?b2cos2tdt，计算卫星轨道周长，并计算卫星运动

的速度与第一宇宙速度（7.9km/秒）比较。如果取地球长半轴为R=6378 （km），卫星速度会如何变化？与第一宇宙速度相比如何？

九 线性规划

例1：某工厂生产A，B两种产品，所用原料均为甲、乙、丙三种：生产一件产品所需原料和所获利润以及库存原料情况如下所示： 原料甲（公斤） 原料乙（公斤） 原料丙（公斤） 利润（元） 产品A 产品B 8 6 4 8 4 6 7000 10000 380 300 220 库存原料量 在该厂只有表中所列库存原料的情况下，如何安排A，B两种产品的生产数量可以获得最大利润？

设生产A产品x1件，生产B产品x2件，z为所获利润，我们将问题归结为如下的线性规划问题：

目标函数：min ?z??(7000x1?10000x2),

（为求z最大，则?z最小，MATLAB默认求最小）

?8x1?6x2?380?约束条件为：?4x1?8x2?300，

?4x?6x?2202?1求解：为方便解释MATLAB求线性规划问题的命令，先将线性规划一般标准形式陈述如下：