注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

倚楼听风雨

没有理想的人,永远也不能翱翔与蓝天白云之上~

 
 
 

日志

 
 

数值分析(matlab)  

2007-11-11 00:13:35|  分类: Mtalab |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |

网上转载
数值分析(1)

 

每当难以对一个函数进行积分、微分或者解析上确定一些特殊的值时,就可以借助计算机在数值上近似所需的结果。这在计算机科学和数学领域,称之为数值分析。至此,可以猜到,MATLAB提供了解决这些问题的工具。本章将介绍这些工具的使用。

13.1  绘图

说到绘图,只要计算函数在某一区间的值,并且画出结果向量,这样就得到了函数的图形。在大多数情况下,这就足够了。然而,有时一个函数在某一区间是平坦的并且无激励,而在其它区间却失控。在这种情况下,运用传统的绘图方法会导致图形与函数真正的特性相去甚远。MATLAB提供了一个称为fplot的巧妙的绘图函数。该函数细致地计算要绘图的函数,并且确保在输出的图形中表示出所有的奇异点。该函数的输入需要知道以字符串表示的被画函数的名称以及2元素数组表示的绘图区间。例如:

>>fplot( humps , [0  2])

>>title( FPLOT OF HUMPS )

02之间计算函数humps,并显示该函数的图形。(见图13.1)

13.1  函数humps的图形

在这个例子中, humps MATLABM文件函数。

function y=humps(x)

%  HUMPS A function used by QUADDEMO, ZERODEMO and FPLOTDEMO.

%  HUMPS(X) is a function with strong maxima near x= .3 and x= .9.

%  See QUADDEMO,  ZERODEMO and FPLOTDEMO.

%  Copyright (c) 1984-93 by The MathWorks, Inc.

y=1 ./ ((x - .3) .^ 2+ .01)+1 ./ ((x - .9) .^ 2+ .04) - 6;

fplot适用于任何具有单输入和单输出向量的函数M文件。即如同humps,输出变量y返回一个与输入x同样大小的数组,在数组到数组意义上yx有联系。在使用fplot(以及其它数值分析函数)的过程中,最普遍犯的错误是忘记把函数名加上引号。即fplot需要知道字符串形式的函数名。如果输入fplot(humps , [0 , 2])MATLAB认为humps是工作空间中的一个变量,而不是函数的名称。注意把变量humps定义为所需要的字符串,就可避免这个问题。

>>humps= humps ;

>>fplot(hump , [0  2])

这时,MATLAB从变量humps中获得字符串 humps

对于可表示成一个字符串的简单的函数,如 fplot绘制这类函数的曲线时,不用建立M文件,只需把x当作自变量,把被绘图的函数写成一个完整的字符串。

>>f= 2*exp(-x) .* sin(x) ;

式中,运用数组乘法定义了函数

>>fplot(f , [0  8]);

>>title(f) , xlabel(x)

13.2 的曲线

在区间 绘出上述函数,产生如图13.2所示的图形。

除了这些基本特性,函数fplot还有很多强大的功能,有关详细的信息,参阅《MATLAB参考指南》或在线帮助。

13.2  极小化

作图除了提供视觉信息外,还常常需要确定一个函数的其它更多的特殊属性。在许多应用中,特别感兴趣的是确定函数的极值,即最大值(峰值)和最小值(谷值)。数学上,可通过确定函数导数(斜率)为零的点,解析上求出这些极值点。检验humps的图形在峰值和谷值点上的斜率就很容易理解这个事实。显然,如果定义的函数简单,则这种方法常常奏效。然而,即使很多容易求导的函数,也常常很难找到导数为零的点。在这种情况下,以及很难或不可能解析上求得导数的情况下,必须数值上寻找函数的极值点。MATLAB提供了两个完成此功能的函数fminfmins。这两个函数分别寻找一维或n维函数的最小值。这里仅讨论fmin。有关fmins的详细信息,参阅《MATLAB参考指南》。因为f(x)的最大值等于-f(x)的最小值,所以,上述fminfmins可用来求最大值和最小值。如果还不清楚,把上述图形倒过来看,在这个状态下,峰值变成了谷值,而谷值则变成了峰值。

为了解释求解一维函数的最小值和最大值,再考虑上述例子。从图13.2可知,在xmax=0.7附近有一个最大值,并且在xmin=4附近有一个最小值。而这些点的解析值为: 。为了方便,用文本编辑器编写一个脚本M文件,并用fmin寻出数值上极值点,给出函数主体如下:

%  ex_fmin.m

fn= 2*exp(-x)*sin(x) ;                 %  define function for min

xmin=fmin(fn , 2 , 5)                    %  search over range 2<x<5

emin=5*pi / 4-xmin                      %  find error

x=xmin;                                      %  need x since fn has x as its variable

ymin=eval(fn)                             %  evaluate at xmin

fx=-2*exp(-x)*sin(x);               %  define for max:note minus sign

xmax=fmin(fx , 0 , 3)                   %  search over range 0<x<3

emax=pi / 4-xmax                        %  find error

x=xmax;                                      %  need x since fn has x as its variable

ymax=eval(fn)                             %  evaluate at xmax

下面是M文件的运行结果:

>>ex-fmin

xmin =

    3.9270

emin =

  1.4523e-006

ymin =

   -0.0279

xmax =

    0.7854

emax =

 -1.3781e-005

ymax =

    0.6448

这些结果与上述图形非常吻合。注意,fmin的工作方式很像fplot。要计算的函数可用一个函数M文件表达,或者只给出一个x为自变量的字符串。上述例子就是使用后一种方法。这个例子也使用了函数eval,它获取一个字符串,并解释它,如同在MATLAB提示符下输入该字符串。由于要计算的函数以x为自变量的字符串形式给出,那么设置x等于xminxmax,允许eval计算该函数,找到yminymax

最后,特别注意,求数值上的最小值包含一个搜索过程,fmin不断计算函数值,寻求其最小值。如果计算的函数需要很大的计算量,或者该函数在搜索区间不止一个最小值,则该搜索过程所花的时间比较长。在有些情况下,搜索过程根本找不到结果。当fmin找不到最小值时,它会停止运行并提供解释。

与函数fmin一样,函数fmins搜索最小值。不过,fmins搜索向量的标量函数的最小值。即fmins寻找

这里x是函数f(.)的向量参数,函数f(.)返回标量值。函数fmins利用单纯形法求最小值,它不需要精确的梯度计算。任何一种优化工具箱中具有更多扩展的优化算法

13.3  求零点

正如人们对寻找函数的极点感兴趣一样,有时寻找函数过零或等于其它常数的点也非常重要。一般试图用解析的方法寻找这类点非常困难,而且很多时候是不可能的。在上述函数humps的图中(如图13.3所示),该函数在x=1.2附近过零。

13.3 humps函数的图形

MATLAB再一次提供了该问题的数值解法。函数fzero寻找一维函数的零点。为了说明该函数的使用,让我们再运用humps例子。

>>xzero=fzero( humps , 1.2)      %  look for a zero near 1.2

xzero=

    1.2995

>>yzero=humps(xzero , 1.2)         %  evaluate at xzero

yzero=

   3.5527e-15

所以,humps的零点接近于1.3。如前所述,寻找零点的过程可能失败。如果fzero没有找到零点,它将停止运行并提供解释。

当调用函数fzero时,必须给出该函数的名称。但由于某种原因,它不能接受以x为自变量的字符串来描述的函数。这样,即使在fplotfmin中都具有的这个特性,fzero将不工作。

fzero不仅能寻找零点,它还可以寻找函数等于任何常数值的点。仅仅要求一个简单的再定义。例如,为了寻找f(x)=c的点,定义函数g(x)=f(x)-c,然后,在fzero中使用g(x),就会找出g(x)为零的x值,它发生在f(x)=c时。

13.4  积分

一个函数的积分或面积也是它的另一个有用的属性。MATLAT提供了在有限区间内,数值计算某函数下的面积的三种函数:trap2 , quadquad8。函数trapz通过计算若干梯形面积的和来近似某函数的积分,这些梯形如图13.4所示,是通过使用函数humps的数据点形成。

13.4 粗略的梯形逼近曲线下的面积示意图

从图中可明显地看出,单个梯形的面积在某一段欠估计了函数真正的面积,而在其它段又过估计了函数的真正面积。如同线性插值,当梯形数目越多时,函数的近似面积越准确。例如,在图13.4中,如果我们大致增加一倍数目的梯形,我们得到如下页(如图13.5)所示的更好的近似结果。

13.5 较好的梯形逼近曲线下的面积示意图

对如上所示的两个曲线,用trapz在区间-1<x<2上计算y=humps(x)下面的面积:

>>x=-1 : 0.17 : 2;                             %  rough approximation

>>y=humps(x);

>>area=trapz(x , y)                          %  call trapz just like the plot command

area =

   25.9174

>>x=-1 : 0.07 : 2;                             %  better approximation

>>y=humps(x);

>>area=trapz(x , y)

area =

   26.6243

自然地,上述两个结果不同。基于对图形的观察,粗略近似可能低估了实际面积。除非特别精确,没有准则说明哪种近似效果更好。很明显,如果人们能够以某种方式改变单个梯形的宽度,以适应函数的特性,即当函数变化快时,使得梯形的宽度变窄,这样就能够得到更精确的结果。

MATLAB的函数quadquad8是基于数学上的正方形概念来计算函数的面积,这些积分函数的操作方式一样。为获得更准确的结果,两个函数在所需的区间都要计算被积函数。此外,与简单的梯形比较,这两个函数进行更高阶的近似,而且quad8quad更精确。这两个函数的调用方法与fzero相同,即

>>area=quad( humps , -1 , 2)  %  find area between -1 and 2

area =

   26.3450

>>area=quad8( humps , -1 , 2)

area =

   26.3450

注意,这两个函数返回完全相同的估计面积,而且这个估计值在两个trapz面积的估计值之间。有关MATLAB的积分函数的其它信息,参阅《MATLAB参考指南》或在线帮助。

13.5  微分

与积分相反,数值微分非常困难。积分描述了一个函数的整体或宏观性质,而微分则描述一个函数在一点处的斜率,这是函数的微观性质。因此积分对函数的形状在小范围内的改变不敏感。而微分却很敏感。一个函数小的变化,容易产生相邻点的斜率的大的改变。

由于微分这个固有的困难,所以尽可能避免数值微分,特别是对实验获得的数据进行微分。在这种情况下,最好用最小二乘曲线拟合这种数据,然后对所得到的多项式进行微分。或用另一种方法,对该数据进行三次样条拟合,然后寻找如第11章所讨论的样条微分。例如,再次考虑第11章曲线拟合的例子。

>>x=[0  .1  .2  .3  .4  .5  .6  .7  .8  .9  1]

>>y=[-.447  1.978  3.28  6.16  7.08  7.34  7.66  9.56  9.48  9.30  11.2];  % data

>>n=2;  %  order of fit

>>p=polyfit(x , y , n)   %  find polynomial coefficients

p =

   -9.8108   20.1293   -0.0317

>>xi=linspace(0 , 1 , 100);

>>z=polyval(p , xi);    %  evaluate polynomial

>>plot(x , y , o ' , x , y , xi , z , ' : ')

>>xlabel( x ) , ylabel( y=f(x) ) , title(Second Order Curve Fitting )

在这种情况下,运用多项式微分函数polyder求得微分。

>>pd=polyder(p)

pd =

  -19.6217   20.1293

13.6 二次曲线拟合

<[if !supportEmptyParas]><[endif]> 

<![endif]> 的微分是dy/dx=-19.6217x+20.1293。由于一个多项式的微分是另一个低一阶的多项式,所以还可以计算并画出该函数的微分。

>>z=polyval(pd , xi);  %  evaluate derivative

>>plot(xi , z)

>>xlabel( x ) , ylabel( dy/dx ) , title(Derivative of a curve Fit Polynimial )

(微分曲线如图13.7所示)

13.7 曲线拟合多项式微分

<[if !supportEmptyParas]><[endif]> 

在这种情况下,拟合的多项式为二阶,使其微分为一阶多项式。这样,微分为一条直线,它意味该微分与x成线性变化。

给定一些描述某函数的数据,MATLAB提供了一个计算其非常粗略的微分的函数。这个函数命名为diff,它计算数组中元素间的差分。因为微分定义为:

y=f(x)的微分可近似为:

        这里h>0

它是y的有限差分除以x的有限差分。因为diff计算数组元素间的差分,所以在MATLAB中,可近似求得函数的微分。继续前一个例子:

>>dy=diff(y) ./ diff(x);      %  compute differences and use array division

>>xd=x(1 : length(x)-1);     %  create new x axis since dy is shorter than y

>>plot(xd , dy);

>>title( Approximate Derivative Using DIFF )

>>ylabel( dy/dx ) , xlabel( x )

13.8 diff得到的近似微分

由于diff计算数组元素间的差分,所以,其所得输出比原数组少了一个元素。这样,画微分曲线时,必须舍弃x数组中的一个元素。当舍弃x的第一个元素时,上述过程给出向后差分近似,而舍弃x的最后一个元素,则给出向前差分近似。比较上述两条曲线,显而易见,用有限差分近似微分会导致很差的结果,特别是被噪声污染了的数据。

  评论这张
 
阅读(1073)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017