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

MATLAB数学实验第四章函数和方程第四章函数和方程4.1预备知识：零点、极值和最小二乘法4.2函数零点、极值和最小二乘拟合的MATLAB指令4.3计算实验：迭代法4.4建模实验：购房贷款的利率4.1预备知识：零点非线性方程f(x)=0•若对于数有f()=0,则称为方程的解或根，也称为函数f(x)的零点•若f()=0,f’()0则称为单根。•若有k1,f()=f’()=…=f(k-1)()=0，但f(k)()0,称为k重根•非线性方程求解通常用数值方法求近似解非线性方程（组）f(x)=0，x=(x1,x2,…,xn),f=(f1,f2,…,fm)4.1预备知识：极值如果对于包含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上的一个全局极小(大)值点。设x为标量或向量，y=f(x)是xD上的标量值函数。4.1预备知识：极值4.1预备知识：最小二乘拟合成绩选手国籍日期13秒24米尔布恩美国1972年9月7日13秒21卡萨那斯古巴1977年8月21日13秒16内赫米亚美国1979年4月14日13秒00内赫米亚美国1979年5月6日12秒93内赫米亚美国1981年8月19日12秒92金多姆美国1989年8月16日12秒91杰克逊英国1993年8月20日12秒88刘翔中国2006年7月12日12秒87罗伯斯古巴2008年6月12日数据：男子110米栏记录19701975198019851990199520002005201012.8512.912.951313.0513.113.1513.213.2513.3问题：估计什么时候突破12秒50？4.1预备知识：最小二乘拟合4.1预备知识：最小二乘拟合•%考虑负指数函数，Matlab程序clear;closeall;data=[13.2413.2113.1613.0012.9312.9212.9112.8812.87];year=[197219771979197919811989199320062008];plot(year,data,'o');y=data;x=year-1972;fun=@(c,x)c(1)*exp(-c(2)*x);%拟合方程y=a*exp(b*x)c=lsqcurvefit(fun,[13,(13.24-12.87)/36],x,y)%结果c(1)=13.1539,c(2)=0.0007fun2=@(x)c(1)*exp(-c(2)*x)-12.5;x2=fsolve(fun2,50)%结果x2=71.4564，约71年以后，即2043年xx=0:75;yy=fun(c,xx);figure;plot(year,data,'o',xx+1972,yy);gridon;4.1预备知识：最小二乘拟合19701980199020002010202020302040205012.412.512.612.712.812.91313.113.213.313.44.1预备知识：最小二乘拟合•假设已知经验公式y=f(c,x)(c为参数,x为自变量),要求根据一批有误差的数据(xi,yi),i=0,1,…,n,确定参数c.这样的问题称为数据拟合。•最小二乘法就是求c使得均方误差最小化Q(c)=•当f关于c是线性函数,问题转化为一个线性方程组求解。•如果f关于c是非线性函数，问题转化为函数极值问题niiixcfy02)),((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中一个多项式用系数降幂排列向量来表示。例2.用2次多项式拟合下列数据.x0.10.20.150-0.20.3y0.950.840.861.061.500.72M文件eg4_2.m例1.求多项式x3+2x2-5的根p=[120-5];x=roots(p),polyval(p,x)返回•x=[0.1,0.2,0.150-0.20.3];•y=[0.950.840.861.061.500.72];•p=polyfit(x,y,2)p=1.7432-1.69591.0850•fun=poly2sym(p,x)%把多项式p转成符号表达式fun=3387/1943*x^2-7637646031980105/4503599627370496*x+4886217849135065/4503599627370496•vpa(fun,5)ans=1.7432*x^2-1.6959*x+1.0850•xi=-0.2:0.01:0.3;•yi=polyval(p,xi);•plot(x,y,‘ro’,xi,yi)%拟合效果作图更简单的方法可以用曲线拟合工具箱cftoolFun=inline(‘funstr’,’var’)定义一个Inline函数,其中funstr是函数的表达式,var是变量名Fun=@Mfun定义一个函数句柄，这里Mfun是函数的M文件表达方式Fun=@(var)funstr定义匿名函数，其中var是变量名,funstr是函数的表达式4.2非线性函数的MATLAB表达x=fzero(Fun,x0)返回一元函数Fun的一个零点,其中Fun为函数句柄、inline函数或匿名函数。x0为标量时,返回函数在x0附近的零点；x0为向量[a,b]时,返回在[a,b]中的零点4.2函数零点MATLAB指令[x,f,h]=fsolve(Fun,x0)x:返回多元函数Fun在x0附近的一个零点，其中x0为迭代初值向量;f:返回Fun在x的函数值,应该接近0;h:返回值如果大于0，说明计算结果可靠，否则计算结果不可靠。例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)例4求方程组在原点附近的解2141101408xxyexyxxx(1)yx(2)functionf=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,[00])使用函数句柄•使用Inline函数fun=inline('[4*x(1)-x(2)+exp(x(1))/10-1,-x(1)+4*x(2)+x(1)^2/8]','x');[x,f,h]=fsolve(fun,[00])•使用匿名函数fun=@(x)[4*x(1)-x(2)+exp(x(1))/10-1,-x(1)+4*x(2)+x(1)^2/8];[x,f,h]=fsolve(fun,[00])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均为向量。4.2函数极值MATLAB指令例5.(1)求函数y=x2+2x+1在[-11]上的极小值点及极小值.(2)当a=2时，求函数y=(x+a)e|x+a|在[-a-1a+1]上的极小值.(3)求二元函数f(x,y)=100(y-x2)2-(1-x)2的极小值。(1)fun1=inline('x.^2+2*x+1','x')fun=Inlinefunction:fun1(x)=x.^2+2*x+1[x,y]=fminbnd(fun1,-1,1)x=-1y=0(2)clearfun2=inline('(x+a).*exp(abs(x+a))','x','a')fun=Inlinefunction:fun2(x,a)=(x+a).*exp(abs(x+a))a=2;[x,y]=fminbnd(@(x)fun2(x,a),-a-1,a+1)%注意含参数形式的调用格式x=-3y=-2.7183•(3)•fun3=inline('100*(v(2)-v(1)^2)^2+(1-v(1))^2','v')fun3=Inlinefunction:fun3(v)=100*(v(2)-v(1)^2)^2+(1-v(1))^2•[v,fv]=fminsearch(fun3,[11])v=11fv=0注:(1)在上述命令中，fun可以用内联函数，匿名函数或编写M函数文件。(2)在使用fsolve,fminsearch等指令时，多变量必须合写成一个向量变量，如用v(1),v(2),…。(3)如何求解函数的极大值？(4)使用帮助文件学习fminunc的使用.4.2最小二乘拟合MATLAB指令•假设已知经验公式y=f(c,x)(c为参数,x为自变量),要求根据一批有误差的数据(xi,yi),i=0,1,…,n,确定参数c.这样的问题称为数据拟合。•最小二乘法就是求c使得均方误差最小化Q(c)=•当f关于c是线性函数,问题转化为一个线性方程组求解，且其解存在唯一。•如果f关于c是非线性函数，问题转化为函数极值问题niiixcfy02)),((c=lsqnonlin(Fun,c0)non-linearleastsquaresproblems.使用迭代法搜索最优参数c.其中Fun是以参数c(可以是向量)为自变量的函数,表示误差向量y-f(c,x)(x,y为数据)，c0为参数c的近似初值(与c同维向量),具体使请看帮助文件。c=lsqcurvefit(Fun2,c0,x,y)从外部输入数据，这里Fun2为两变量c和x的函数f(c,x)用lsqnonlin解例4.2•先写M函数fitf.mfunctione=fitf(c)x=[0.10.20.150-0.20.3];y=[0.950.840.861.061.500.72];e=y-(c(1)*x.^2+c(2)*x+c(3));•然后执行c=lsqnonlin(@fitf,[0,0,0])%这里[0,0,0]为c的预估值，作为迭代初值。例4.2用lsqcurvefit解例4.2fun2=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.69581.0850•注意lsqnonlin与lsqcurvefit用法上的区别迭代法是从解的初始近似值x0(简称初值)开始，利用某种迭代格式xk+1=g(xk),求得一近似值序列x1,x2,…,xk,xk+1,…逐步逼近于所求的解(称为不动点)。最常用的迭代法是牛顿迭代法，其迭代格式为xxfxfxkkkk1()'()1迭代法4.3计算实验：迭代法牛顿法程