《应用实例仿真与高效算法开发》由会员分享，可在线阅读，更多相关《应用实例仿真与高效算法开发（39页珍藏版）》请在金锄头文库上搜索。
1、第三章 微积分问题的数值实验通过本章的介绍， 加上一些具有代表性的例题， 读者可以了解微积分学的一些重要的概 念与方法。实际科学与工程研究中，往往只得到一些离散的实验数据，并不知晓函数本身， 无法用解析的方法对这些数据进行处理， 此时，通过数值的方式进行数值微分与数值积分的 运算显得尤为重要。本章的内容是围绕着单变量与多变量函数微积分、 函数极限、 级数求和、 Taylor 幂级数 展开 Fourier 级数展开等问题展开的。为方便读者理解，每节中都附有代表性的例题，并给 出了程序代码， 读者可以根据自己的需要对程序进行扩充， 增补，以实现自己所需要的功能。 作者的主要用意是综合运用不同指令解
2、决具体问题， 为读者解决实际问题提供一些思路和借 鉴。本章的主要内容包括微积分问题的解析解函数的级数展开与求和数值微分问题数值积分问题曲线积分与曲面积分的计算3.1微积分问题的解析解3.1.1极限问题的解析解函数极限问题的定义为：设函数f(x)在点X。的某一去心领域内有定义，如果对于任意给定的正数(无论它多么小)，总存在正数，使得对于适合不等式Ov X-X。的一切X，对应的函数值f (x)都满足不等式f(x)-A X0时的极限，记作XimXL f(X)= A 或 f(X)T A (当 XT X0)上述定义中的Xo可以是某确定的值，也可为无穷大。当常数A满足：A=f(x0)时，即密 f(x)=f
3、(x0)时，可以称函数f (x)在点xo连续。函数在某一点连续又可分为左连续和右连续。如果lim f (x)二 f (xo - 0)JXo -0存在且等于f(Xo)，就说函数f (X)在点处左连续。如果xlirX0 Of(x f(X00)存在且等于f(Xo)，就说函数f (x)在点处右连续。极限问题在MATLAB符号运算工具箱中提供了很方便的指令limit()，该通过该函数不同的调用格式，可分别求得函数的极限和左/右极限。其调用格式如下：P = limit(fun, x, xo)求函数fun关于自变量x在xo处的极限P = limit(fun, x, xo， left或right) 求函数fu
4、n关于自变量x在Xo处的单边极限下面通过例子来演示MATLAB求函数极限的方法。2【例3-1】求解极限问题ln(a bx )secx 一 cosx分析：Matlab中求解此类问题首先要定义符号变量，然后定义极限式子，接着才调用 函数球给定函数的极限。 syms a b x; fun=log(a+b*xA2)/(sec(x)-cos(x); Lim = limit(fun, x, 1)Lim =-Iog(a+b)*cos/(-1+cos A2)【例3-2】求解单边极限问题lim xsin(x)。 X严 syms x; fun=xA(sin(x); Iim=Iimit(fun,x,0,right)
5、lim =1我们还可以绘制出(-0.1,0.1)区间的函数曲线，如下图所示 x=-0.1:0.0001:0.1; f_x=x.A(sin(x);plot(x,f_x,black-,0,1,blacko)图3-1x=0附近的曲线通过看x关于原点邻域内的曲线可以看到，函数f (x)二x(sinx)在x=0处是连续的，所以无论是求左极限还是右极限，都为1。另外，上面我们接触的都是函数自变量为单个情形，接下来简单介绍一下多变量函数的极限问题。多元函数的极限问题同样可以通过嵌套调用matlab指令limit()来解决。如面对二元函数f(x,y)，若想求得二元函数的极限Lim 二 lim f (x, y)
6、xoy 砂o嵌套调用limit()函数：lim = limit(limit( f ,x,x0), y, y0)也可以这样嵌套调用:lim =limit(limit( f, y, y),x, x0x_；(y)，则极限【例3-3】求二元函数的极限问题limx-：y2-3/(x2+y2)sin Xe2X(12ax-by值得指出的是，当，不是确定的值，而是另外一个自变量的函数，例如 的求取顺便不能交换。 syms a b x y fun=exp(-3/(xA2+yA2)*sin(x)A2/xA2*(1+1/yA2)A(a*x+bA2*yA2) Lim=limit(limit(fun,x,1/sqrt(
7、y),y,inf)Lim =exp(bA2)此例题中如x 某一确定数，那么极限问题中X ,y的顺序可以调换。3.1.2函数导数的解析解数学上可以表达导数问题的物理意义可以理解为非匀速直线运动的速度和切线的斜率,为：limX %f (x) - f (Xo)X Xo更严格的定义为，设函数 y二f(x)在某个邻域内有定义，当自变量x在xo处取得增量 时，相应地函数y取得增量 y = f (x x f (x0)；如果与之比当=x 0时的极限存在, 则称函数y = f (x)在点xo处可导，并称这个极限为函数 y二f (x)在点xo处的导数，记为 y ，即X丸?耳.3.f(X。+3) f(X。)ylim
8、 lim -XzXo.Xx.X -xMATLAB中提供的求函数导数的指令为diff()，它可以解出给定函数的各阶导数，其调用格式为：y = diff (fun ,x)求函数fun的关于x的导数y = diff (fun, x,n)求函数fun的关于x的 n阶导数说明：其中fun为给定的待求导函数，x为自变量，这两个变量都为符号型的，求导的 阶次用n来表示，默认为一阶导数。下面举例说明。【例3-4】 给定函数f(x)二3COSX，试求其关于自变量 x的四阶导数。x3+7x+2首先看看一阶导数 syms x; f=cos(x)/(xA3+7*x+2); f_1= diff(f,x)2sin(x)c
9、os(x) (3 x + 7)332x + 7 x + 2(x + 7 x + 2) pretty(f_1)这样以来增加了结果的可读性。一阶导数在每个点上值与原函数的值可以通过 MATLAB函数绘制出来如下图所示：x_1=0:0.001:5;y=subs(f,x,x_1); y_1=subs(f_1,x,x_1);plot(x_1,y,x_1,y_1,:)图3-2函数及其一阶导数 F面求原函数的四阶导数。 f_4=diff(f,x,4); latex(f4) latex(f_4)得出的结果冗长，用latex()命令可以得到更好的显示效果：cos(x)sin(x)(3x2 7)cos(x)(3x
10、2 7)2cos(x)xsin (x)(3x2 7)3x3+7x+2(x3 7x 2)2(x3 - 7x - 2)3(x3 7x - 2)2 (x3 - 7x 2)4144sin (x)(3x27)x(x3 7x 2)224sin (x)(x3 7x 3)224cos(x)(3x27)4(x3 7x 2)5 216cos(x)(3x27)2x(x3 7x 2)4216cos(x)x(x3 7x 2)348cos(x)(3x27)(x3 7x 2)3为达到更好的显示效果，我们选择了另外的化简方法：分别执行collect(simple(f_4),cos(x) 和 collect(simple(f_
11、4),sin(x)。此外，matlab中提供的diff()指令可以对多元函数求偏导数。调用格式如下，假设已知二元函数f(x, y)，欲求fC：xyn)，则可以通过以下函数求得，f = diff (diff (f, x, m), y, n)，或者 f = diff (diff (f, y, n), x, m)22【例3-5】求二元函数fun = f（X,y） =（y2-3y）ex3xy）的偏导数，并用图形表示出来。 syms x y; fun=(yA2-3*y)*exp(-(xA2+yA2+3*x*y); fun_x=simple(diff(fun,x) fun_x=diff(fun,x)fun
12、_x =(yA2-3*y)*(-2*x-3*y)*exp(-xA2-yA2-3*x*y) fun_y=diff(fun,y)fun_y =(2*y-3)*exp(-xA2-yA2-3*x*y)+(yA2-3*y)*(-2*y-3*x)*exp(-xA2-yA2-3*x*y) x,y=meshgrid(-3:0.1:3,-2:0.1:2); fun=(y.A2-3*y).*exp(-(x.A2+y.A2+3*x.*y);fun_x=(y.A2-3*y).*(-2*x-3*y).*exp(-x.A2-y.A2-3*x.*y);fun_y=(2*y-3).*exp(-x.A2-y.A2-3*x.*y
13、)+(y.A2-3*y).*(-2*y-3*x).*exp(-x.A2-y.A2-3*x.*y); figure(1); surf(x,y,fun)1500 图3-3原二元函数的三维图形接着察看对X的偏导数， figure(2); surf(x,y,fun_x)再看对y的偏导数， figure(3); surf(x,y,fun_y)IDCC-图3-4原函数对x的偏导数图3-5原函数对y的偏导数对隐函数和参数方程的导数问题同样可通过diff()指令来求取。举例，【例 3-6】对二元函数 fun 二 f (x,y) = (y? -3y)ex 审 3xy)，求。经分析可知，訶;:x_fun(x, y
14、)fun (x, y)，其中分子和分母上的偏导数可以通过x/cydiff()很方便的求得。 syms x y; fun=(yA2-3*y)*exp(-(xA2+yA2+3*x*y);Y_x = simple(diff(fun,x)/diff(fun,y)Y_x =-y*(y-3)*(2*x+3*y)/(-2*y+3+2*yA3+3*yA2*x-6*yA2-9*x*y)为便于阅读，表达式可写为：-y(y - 3)(2x +3y)(2y3 3xy2 - 6y2 - 2y -9xy)对于参数方程y二f(d),x二g(d)那么，y关于x的k阶偏导数可调用 diff()命令求得,具体计算格式如下：dkyd y 二diff (fj,k);diff (gjk) dx22【例 3-7】y = -3te4,x 二 e，试求，d2y dx2。 syms t y=-3*t*exp(-tA2-1); x=exp(-tA2-1);simple(diff(y,t,2)/diff(x,t,2)ans =(9*t-6*tA3
