MATLAB数学实验 第四章 函数和方程

第四章 函数和方程 4.1 预备知识：零点、极值和最小二乘法 4.2 函数零点、极值和最小二乘拟合的MATLAB指令 4.3 计算实验：迭代法 4.4 建模实验：购房贷款的利率

f (x) = 0， x=(x1, x2, …, xn), f=(f1, f2, …, fm) 4.1 预备知识：零点 非线性方程 f (x) = 0 若对于数有f () = 0, 则称为方程的解或根，也称为函数f (x)的零点 若f () = 0, f ’()0 则称为单根。 若有k >1, f () = f ’() = …= f (k-1)() = 0，但 f (k)()0 , 称为k重根 非线性方程求解通常用数值方法求近似解 非线性方程（组） f (x) = 0， x=(x1, x2, …, xn), f=(f1, f2, …, fm)

4.1 预备知识：极值 设x为标量或向量，y=f(x)是xD上的标量值函数。 如果对于包含x=a的某个邻域 ,有 f(a)f(x) (f(a)f(x))对任意x成立, 则称a为f(x)的一个局部极小(大)值点。 如果对任意xD，有f(a)f(x)（f(a)f(x)）成立， 则称a为f(x)在区域D上的一个全局极小(大)值点。

4.1 预备知识：极值

4.1 预备知识：最小二乘拟合 假设已知经验公式y=f(c,x)(c为参数, x为自变量), 要求根据一批有误差的数据(xi,yi), i=0,1,…,n, 确定参数c.这样的问题称为数据拟合。 最小二乘法就是求c使得均方误差最小化 当f关于c是线性函数,问题转化为一个线性方程组求解，且其解存在唯一。 如果f关于c是非线性函数，问题转化为函数极值问题

4.2 函数零点MATLAB指令 多项式 y=polyval(p,x) 求得多项式p在x处的值y，x可以是一个或多个点 p3=conv(p1,p2) 返回多项式p1和p2的乘积 [p3,r]=deconv(p1,p2) p3返回多项式p1除以p2的商，r返回余项 x=roots(p) 求得多项式p的所有复根. p=polyfit(x,y,k)用k次多项式拟合向量数据(x, y),返回多项式的降幂系数

MATLAB中一个多项式用系数降幂排列向量来表示。 例1.求多项式x3 + 2 x2 - 5的根 » p=[1 2 0 -5]; x=roots(p) , polyval(p,x) 例2.用2次多项式拟合下列数据. x 0.1 0.2 0.15 0 -0.2 0.3 y 0.95 0.84 0.86 1.06 1.50 0.72 M文件Ch4_2_1ex2.m

4.2 非线性函数的MATLAB表达 Fun=inline(‘funstr’,’var’) 定义一个inline Fun=@Mfun 定义一个函数句柄，这里Mfun是函数的M文件表达方式 Fun=inline(‘funstr’,’var’) 定义一个inline 函数,其中funstr是函数的表达式, var是变量名 Fun=@(var)funstr 定义匿名函数，其中var是变量名, funstr是函数的表达式

4.2 函数零点MATLAB指令 x=fzero(Fun, x0) 返回一元函数Fun的一 个零点,其中Fun为函数句柄、 inline函数或匿名函数。 x0为标量时, 返回函数在x0附近的零点； x0为向量[a, b]时, 返回在[a,b]中的零点

[x,f,h]=fsolve(Fun, x0) x: 返回多元函数Fun在x0附近的一个零点，其中x0为迭代初值向量; f: 返回Fun在x的函数值, 应该接近0; h: 返回值如果大于0，说明计算结果可靠，否则计算结果不可靠。

>> [x,f,h]=fsolve(fun,-0.6) 例3 求函数y=xsin(x2-x-1)在(-2, -0.1)内的 零点 >> fun=inline('x*sin(x^2-x-1)','x') >> fzero(fun,[-2 -0.1]) >> fzero(fun,[-2,-1.2]), fzero(fun,[-1.2,-0.1]) >> fzero(fun,-1.6), fzero(fun,-0.6) >> [x,f,h]=fsolve(fun,-1.6) >> [x,f,h]=fsolve(fun,-0.6)

f(1)=4*x(1)-x(2)+exp(x(1))/10-1; f(2)=-x(1)+4*x(2)+x(1)^2/8; 例4 求方程组在原点附近的解 xx(1) y x(2) function f=eg4_4fun(x) f(1)=4*x(1)-x(2)+exp(x(1))/10-1; f(2)=-x(1)+4*x(2)+x(1)^2/8; >> [x,f,h]=fsolve(@eg4_4fun,[0 0])

4.2 函数极值MATLAB指令 min(y) 返回向量y的最小值 max(y) 返回向量y的最大值 [x,f]=fminbnd(fun,a,b) x返回一元函数y=f(x)在[a,b]内的 局部极小值点, f返回局部极小值 fun为函数句柄或inline函数或匿名函数。 [x,f]=fminsearch(fun,x0) x返回多元函数y=f(x)在初始值x0 附近的局部极小值点,f返回局部极小值. x, x0均为向量。

例 5 .求二元函数f(x,y)= 5-x4-y4+4xy在原点附近的极大值。 解：max fmin(-f) x x(1), y x(2) >> fun=inline(' -5 +x(1)^4+x(2)^4-4*x(1)*x(2)'); >> [x,g]=fminsearch(fun,[0,0]) 注：在使用fsolve, fminsearch等指令时， 多变量必须合写成一个向量变量，如用x(1), x(2),…。

4.2 最小二乘拟合MATLAB指令 假设已知经验公式y=f(c,x)(c为参数, x为自变量), 要求根据一批有误差的数据(xi,yi), i=0,1,…,n, 确定参数c.这样的问题称为数据拟合。 最小二乘法就是求c使得均方误差最小化 当f关于c是线性函数,问题转化为一个线性方程组求解，且其解存在唯一。 如果f关于c是非线性函数，问题转化为函数极值问题

c= lsqnonlin(Fun,c0) 使用迭代法搜索最 优参数c. 其中Fun是以参数c(可以是向量) 为自变量的函数,表示误差向量 y-f(c,x)(x, y为数据)， c0为参数c的近似初值(与c同维向量) c=lsqcurvefit(Fun2,c0, xdata, ydata) 从外部输入数据xdata, ydata， 这里Fun2为两变量c和x的函数 f(c, x)

用lsqnonlin再解例4.2 先写M函数fitf.m function e=fitf(c) x=[0.1 0.2 0.15 0 -0.2 0.3]; y=[0.95 0.84 0.86 1.06 1.50 0.72]; e=y-(c(1)*x.^2+c(2)*x+c(3)); 然后执行 >> c=lsqnonlin(@fitf,[0,0,0]) %这里[0,0,0]为c的预估值，作为迭代初值。

用lsqcurvefit再解例4.2 计算结果 c = 1.7427 -1.6958 1.0850 >>fun2 = inline('c(1)*x.^2+c(2)*x+c(3)','c','x') >> x=[0.1, 0.2, 0.15, 0, -0.2, 0.3]; >> y=[0.95,0.84,0.86,1.06,1.50,0.72]; >> c=lsqcurvefit(fun2,[0,0,0],x,y) 计算结果 c = 1.7427 -1.6958 1.0850 注意lsqnonlin与lsqcurvefit用法上的区别

4.3 计算实验：迭代法 1 迭代法 迭代法是从解的初始近似值x0(简称初值)开始，利用某种迭代格式x k+1 = g (x k ),  1 迭代法 迭代法是从解的初始近似值x0(简称初值)开始，利用某种迭代格式x k+1 = g (x k ), 求得一近似值序列x1, x2, …, xk, xk+1, … 逐步逼近于所求的解(称为不动点)。 最常用的迭代法是牛顿迭代法，其迭代格式为

牛顿法程序newton.m function x=newton(fname,dfname,x0,e) if nargin<4, e=1e-4;end %默认精度x=x0;x0=x+2*e; %使while成立，且进入while后x0得到赋值 while abs(x0-x)>e x0=x; x=x0-feval(fname,x0)/feval(dfname,x0); end

例6 求方程 x 2 - 3 x + e x = 2 的正根 (要求精度 = 10 -6) 解 令f (x) = x 2 - 3 x + e x - 2, f(0)=-1, 当x > 2, f (x) > 0, f ’(x) > 0 即f (x)单调上升，所以根在[0,2]内。 先用图解法找初值， 再用牛顿法程序newton.m求解。 M文件Ch4_3_1ex6.m

2．线性化拟合 例7 .用函数y=aebx 拟合例2的数据 方法一：非线性拟合，记a=c(1), b=c(2) >> fun = inline('c(1)*exp(c(2)*x)','c','x'); >> x=…;y=…; c=lsqcurvefit(fun,[0,0],x,y) 方法二：线性化拟合，两边取对数得 z=lny=lna+bx 转化为线性拟合。 M文件Ch4_3_2ex7.m

4.4 建模实验：购房贷款的利率 例8 .下面是《新民晚报》2000年3月30日上的一则房产广告： 不难算出，你向银行总共借了25.2万， 30年内共要还51.696万， 这个案例中贷款年利率是多少呢？

解 设xk为第k个月的欠款数， a为月还款数, r为月利率。 xk+1 = (1+r) xk- a 那么 xk = (1+r) xk-1- a = (1+r)2 xk-2 – (1+r)a – a =…… = (1+r)k x0 – a[1+(1+r)+……+(1+r)k-1] = (1+r)k x0 – a[(1+r)k-1]/r 根据 a=0.1436, x0=25.2, x360=0 得到 25.2(1+r)360 – 0.1436[(1+r)360-1]/r=0 很难用roots求解！

常识上，r应比当时活期存款月利率略高一些。我们用活期存款月利率0.0198/12 作为迭代初值，用fzero求解 >>clear; fun=inline('25.2*(1+r)^360-((1+r)^360-… 1)/r*0.1436' ,'r') >>r=fzero(fun,0.0198/12); >>R=12*r 得年利率为5.53%.

月还款计算公式 月还款 x0----借款额； N---剩余月数； r---月利率=年利率/12

4.4 建模实验：最佳订货量 每次订货需要收取一定量的生产准备费。 没用完的配件，要在仓库里储存一段时间，为此要付出储存费。 若订货量很小，则需频繁定货，造成生产准备费的增加； 反之，若订货量很大，定货周期延长而使生产准备费减少但会造成储存费的增加。 如何确定合适的订货量？

解 先作一些必要的假设将问题简化 1）汽车工厂对配件的日需求量是恒定的， 每日为r件； 2）所订配件按时一次性交货， 生产准备费每次k1元； 3）储存费按当日实际储存量计算， 储存费每日每件k2元; 4）你的工厂不允许缺货。 设一次订货x件，则订货周期为 T= x/r, 第t天的储存量为 q(t)= x-r t, 0<t<T

第t天的储存费为 k2q(t) 一个周期的总储存费为 一个周期总费用 C(x) = k1+k2x2/(2r) 优化目标是使单位产品费用 f(x)=C(x)/x=k1/x+k2x/(2r) 达到最小 由f’(x)=0 ,即 -k1/x2+k2/(2r)=0 得 （经济批量订货公式）

4.4 建模实验：混沌 线性迭代 xk+1=axk+b 要么收敛于它的不动点x=b/(1-a)，要么趋于无穷大。 不收敛的非线性迭代可能会趋于无穷大，也可能趋于一个周期解， 但也有可能在一个有限区域内杂乱无章地游荡，这类由确定性运动导致的貌似随机的现象称为混沌现象

*平衡与稳定 若g () = ,称为映射g(x)的不动点 若对于不动点附近的初始值x0，迭代收敛于此不动点，称此不动点是稳定的 *昆虫数量的Logistic模型 xk+1 = a x k (1 - x k), 0a4 xk表示第k代昆虫数量(1表示理想资源 环境最大可能昆虫数量)。 a为资源系数 0a4保证了xk在区间(0,1)上封闭。

且由|g’(0)| =a<1,可知它是稳定的, 说明资源匮乏时,昆虫趋于消亡; 对于Logistic模型 解得有两个不动点0和1-1/a 当0a<1, 在[0，1]内有一个不动点0， 且由|g’(0)| =a<1,可知它是稳定的, 说明资源匮乏时,昆虫趋于消亡; 当a>1, 不动点0不再稳定; 当1<a3，由|g’ (1-1/a)| =|2-a|<1可知 不动点1 - a -1稳定，说明资源适当时， 昆虫稳定于一定数量。

a>3，出现两个周期2解，且3<a<1+ 周期2轨道稳定。 迭代开始发生所谓倍周期分岔， 从周期2,周期4,…,周期2 n,… 直到 a = 3.569945672…。 说明a在[3，a]取值时， 昆虫数量呈现规律性振荡。 a >a 迭代序列几乎杂乱无章， 即所谓混沌。

*混沌的特征 (i)初值敏感性: 两个任意近的点出发的两条轨迹 迟早会分得很开; (ii)遍历性: 任意点出发的轨迹总会进入 [0,1]内任意小的开区间。

蛛网图 例9（蛛网图）我们用蛛网图来显示混沌的遍历性。 yk = a x k (1 - x k), xk+1 = yk 蛛网图正好显示迭代计算 x0, y0, x1, y1,……的一系列变化过程。 eg4_9.m

习题 ex1, ex2, ex6, ex8(1) ex9, ex10(1)(3), ex12