应用实例仿真与高效算法开发

1、第三章 微积分问题的数值实验通过本章的介绍，加上一些具有代表性的例题，读者可以了解微积分学的一些重要的概念与方法。实际科学与工程研究中，往往只得到一些离散的实验数据，并不知晓函数本身，无法用解析的方法对这些数据进行处理，此时，通过数值的方式进行数值微分与数值积分的运算显得尤为重要。本章的内容是围绕着单变量与多变量函数微积分、函数极限、级数求和、Taylor幂级数展开Fourier级数展开等问题展开的。为方便读者理解，每节中都附有代表性的例题，并给出了程序代码，读者可以根据自己的需要对程序进行扩充，增补，以实现自己所需要的功能。作者的主要用意是综合运用不同指令解决具体问题，为读者解决实际问题提供

2、一些思路和借鉴。本章的主要内容包括微积分问题的解析解函数的级数展开与求和数值微分问题数值积分问题曲线积分与曲面积分的计算3.1 微积分问题的解析解3.1.1极限问题的解析解函数极限问题的定义为：设函数在点x0的某一去心领域内有定义，如果对于任意给定的正数（无论它多么小），总存在正数，使得对于适合不等式的一切x，对应的函数值都满足不等式那么常数A就叫做函数当时的极限，记作或(当)上述定义中的x0可以是某确定的值，也可为无穷大。当常数A满足：时，即时，可以称函数在点x0连续。函数在某一点连续又可分为左连续和右连续。如果存在且等于，就说函数在点处左连续。如果存在且等于，就说函数在点处右连续。极限问题

3、在MATLAB符号运算工具箱中提供了很方便的指令limit( )，该通过该函数不同的调用格式，可分别求得函数的极限和左/右极限。其调用格式如下：P = limit(fun, x, x0) 求函数fun关于自变量x在x0处的极限P = limit(fun, x, x0，'left' 或 'right') 求函数fun关于自变量x在x0处的单边极限下面通过例子来演示MATLAB求函数极限的方法。【例3-1】求解极限问题。分析：Matlab中求解此类问题首先要定义符号变量，然后定义极限式子，接着才调用函数球给定函数的极限。>> syms a b x;>

4、;> fun=log(a+b*x2)/(sec(x)-cos(x);>> Lim = limit(fun, x, 1)Lim = -log(a+b)*cos(1)/(-1+cos(1)2)【例3-2】求解单边极限问题。>> syms x;>> fun=x(sin(x);>> lim=limit(fun,x,0,'right') lim = 1我们还可以绘制出(-0.1, 0.1)区间的函数曲线，如下图所示>> x=-0.1:0.0001:0.1;>> f_x=x.(sin(x);>>plo

5、t(x,f_x,'black-',0,1,'blacko')图 3-1 x=0 附近的曲线通过看x关于原点邻域内的曲线可以看到，函数在x=0处是连续的，所以无论是求左极限还是右极限，都为1。另外，上面我们接触的都是函数自变量为单个情形，接下来简单介绍一下多变量函数的极限问题。多元函数的极限问题同样可以通过嵌套调用matlab指令limit( )来解决。如面对二元函数，若想求得二元函数的极限嵌套调用limit( )函数：也可以这样嵌套调用：值得指出的是，当，不是确定的值，而是另外一个自变量的函数，例如，则极限的求取顺便不能交换。【例3-3】 求二元函数的极限问题。

6、>> syms a b x y>> fun=exp(-3/(x2+y2)*sin(x)2/x2*(1+1/y2)(a*x+b2*y2)>> Lim=limit(limit(fun,x,1/sqrt(y),y,inf) Lim =exp(b2)此例题中如某一确定数，那么极限问题中x ,y的顺序可以调换。3.1.2函数导数的解析解导数问题的物理意义可以理解为非匀速直线运动的速度和切线的斜率，数学上可以表达为：更严格的定义为，设函数在某个邻域内有定义，当自变量x在x0处取得增量时，相应地函数y取得增量；如果与之比当时的极限存在，则称函数在点x0处可导，并称这个极限

7、为函数在点x0处的导数，记为，即MATLAB中提供的求函数导数的指令为diff( )，它可以解出给定函数的各阶导数，其调用格式为：求函数fun的关于x的导数求函数fun的关于x的n阶导数说明：其中fun为给定的待求导函数，x为自变量，这两个变量都为符号型的，求导的阶次用n来表示，默认为一阶导数。下面举例说明。【例3-4】 给定函数，试求其关于自变量x的四阶导数。首先看看一阶导数>> syms x;>> f=cos(x)/(x3+7*x+2);>> f_1=diff(f,x)>> pretty(f_1) 2 sin(x) cos(x) (3 x +

8、 7) - - - - 3 3 2 x + 7 x + 2 (x + 7 x + 2)这样以来增加了结果的可读性。一阶导数在每个点上值与原函数的值可以通过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_4=diff(f,x,4);>> latex(f4)>> latex(f_4)得出的结

9、果冗长，用latex( )命令可以得到更好的显示效果：为达到更好的显示效果，我们选择了另外的化简方法：分别执行collect(simple(f_4),cos(x)和collect(simple(f_4),sin(x)。此外，matlab中提供的diff( )指令可以对多元函数求偏导数。调用格式如下，假设已知二元函数，欲求，则可以通过以下函数求得，或者【例3-5】 求二元函数的偏导数，并用图形表示出来。>> syms x y;>> fun=(y2-3*y)*exp(-(x2+y2+3*x*y);>> fun_x=simple(diff(fun,x)>&g

10、t; fun_x=diff(fun,x) fun_x = (y2-3*y)*(-2*x-3*y)*exp(-x2-y2-3*x*y) >> fun_y=diff(fun,y) fun_y =(2*y-3)*exp(-x2-y2-3*x*y)+(y2-3*y)*(-2*y-3*x)*exp(-x2-y2-3*x*y)>> x,y=meshgrid(-3:0.1:3,-2:0.1:2);>> fun=(y.2-3*y).*exp(-(x.2+y.2+3*x.*y);>>fun_x=(y.2-3*y).*(-2*x-3*y).*exp(-x.2-y.2

11、-3*x.*y);>>fun_y=(2*y-3).*exp(-x.2-y.2-3*x.*y)+(y.2-3*y).*(-2*y-3*x).*exp(-x.2-y.2-3*x.*y);>> figure(1); surf(x,y,fun)图 3-3 原二元函数的三维图形接着察看对x的偏导数，>> figure(2); surf(x,y,fun_x)再看对y的偏导数，>> figure(3); surf(x,y,fun_y)图 3-4 原函数对x的偏导数 图 3-5 原函数对y的偏导数对隐函数和参数方程的导数问题同样可通过diff( )指令来求取。举

12、例，【例3-6】 对二元函数，求。经分析可知，其中分子和分母上的偏导数可以通过diff( )很方便的求得。>> syms x y;>> fun=(y2-3*y)*exp(-(x2+y2+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*y3+3*y2*x-6*y2-9*x*y)为便于阅读，表达式可写为：对于参数方程那么，y关于x的k阶偏导数可调用diff( )命令求得，具体计算格式如下：【例3-7】，试求，。>> syms t>>

13、; y=-3*t*exp(-t2-1);>> x=exp(-t2-1);>>simple(diff(y,t,2)/diff(x,t,2)ans =(9*t-6*t3)/(-1+2*t2)为便于阅读，可写为：3.1.3积分问题的解析解如果将积分问题分的更细一点的话可分为，不定积分、定积分、无穷积分和多重积分。下面就各类问题分别作简要的说明。.1 不定积分对于可积函数，matlab符号运算工具箱提供的int( )函数可用计算机代替繁重的手工推导，很快捷的得到问题的解。其调用格式如下：该函数得出的为原函数，不定积分的应该是组成的函数族，C为任意常数。考虑函数，用diff()函

14、数求其一阶导数之后，然后对其导数再进行积分，那么通过观察得到的结果就可以验证该命令是否可以得到正确的结果。>> syms x;>> fun=cos(x)/(x2+7*x+6); %定义函数fun>> D_fun=diff(fun); %求fun函数的一阶导数>> R_fun=int(D_fun) %对所得导数进行积分得出原函数 R_fun = 1/5*cos(x)/(x+1)-1/5*cos(x)/(x+6)>>simple(R_fun)ans =cos(x)/(x+1)/(x+6)【例3-8】 试求函数的不定积分。>>

15、syms t %定义符号变量t>> y=(1+t+t2)/t*(1+t2); >> int(y) ans = 1/4*t4+1/3*t3+t2+t+log(t)更直观一点可写为：【例3-9】求函数的不定积分。>> syms t;>> y=exp(t)*(1-exp(-t)/sqrt(t);>> D_y=int(y)D_y= exp(t)-2*t(1/2)>>t=2:0.05:3;>>y1=exp(t)*(1-exp(-t)/sqrt(t);>>D_y1= exp(t)-2*t.(1/2);>&

16、gt;plot(t,y1);>>hold on>>plot(t,D_y1,o)图 3-6 函数y及其定积分曲线此外，对于一些不可积的表达式，matlab可给出一个拟解析解。【例3-10】 求函数的不定积分表达式。>>syms x;>> y=2*exp(-x2/2);>> int(y)ans =pi(1/2)*2(1/2)*erf(1/2*2(1/2)*x)虽然函数y本来不可积，但是matlab中应用数学方法定义了一个符号这样就得到了一个拟解析不定积分表达式。.2 定积分、无穷积分与多重积分对于这三类积分问题，仍然可以用matlab指令

17、int()的调用来解决。定积分和无穷积分可以归结为一类问题，即，当定积分的积分域为无穷大时，便可认为是无穷积分问题。一般的调用格式如下：其中，fun为被积分函数，x为积分变量，m和n分别为积分上下限。除了通用积分指令以外，matlab提供了一类交互式近似积分指令：其中，fun依然为积分函数，可以是字符串或符号形式，a与b为积分区间。它运行后得到一个界面，里面用梯形近似表示积分值，且窗体的右上方的数字是近似积分值，梯形数越大，近似积分的精度逾高。【例3-11】试对，在区间上求定积分。>> syms x>> fun=(x+sin(x)/(1+cos(x);>>

18、int(fun,x,0,pi/2)ans =1/2*pi再用近似积分公式观察积分结果见图3-7，并作一定的比较。>>rsums(fun,0,pi/2)图 3-7 再看一无穷积分的例子，【例3-12】 试对在上求定积分。>> syms x;>> fun=2/exp(x)+3/x5;>> int(fun,x,1,inf) ans = 2*exp(-1)+3/4接下来看一个多重积分的例子，【例3-13】 求积分。>> syms x y z;>> fun=x3+y3+z3;>> Res_int = int(int(in

19、t(fun,z,sqrt(y),x3*y),y,sqrt(x),x3),x,1,2)>>NUM32_Res_int=vpa(Res_int) %积分结果用32位表示NUM32_Res_int =3.2 函数的级数展开与求和多变量或单变量函数函数的级数展开与求和通常包括，Taylor幂级数展开、Fourier级数展开和有/无穷技术求和。本节会针对此类问题的计算机求解方法作一些介绍。 Taylor 幂级数展开对于单变量函数的Taylor幂级数展开，它的数学表达如下：如果在点附近进行Taylor幂级数展开，则有，其中，系数的表达通式为：，如果，那么该幂级数展开又称Maclaurin展开。

20、Taylor幂级数展开可以用符号运算工具箱的taylor( )函数直接导出。其调用格式为 %按进行Taylor幂级数展开其中，fun为函数的符号表达式，x为自变量。该函数将得到项数为k的级数展开表达式。下面举例来演示该指令的应用。【例3-14】 对函数求Talylor幂级数展开式的前7项，并关于和分别对原函数进行Taylor幂级数展开。>> syms x;>> fun=cos(x)/(x2+2*x+1);>> fun_Taylor_expa=taylor(fun,x,8) fun_Taylor_expa = 1-2*x+5/2*x2-3*x3+85/24*x

21、4-49/12*x5+3329/720*x6-1859/360*x7更直观的显示：接下来计算关于的Taylor幂级数展开的前3项，指令如下：>>taylor(fun,x,3,3)ans = 1/16*cos(3)+(-1/16*sin(3)-1/32*cos(3)*(x-3)+(-5/256*cos(3)+1/32*sin(3)*(x-3)2更直观的显示：再观察关于的Taylor幂级数展开，给出下列语句：>> syms h;>> taylor(fun,x,3,h) ans = cos(h)/(h2+1+2*h)+(-sin(h)-cos(h)/(h2+1+2

22、*h)*(2+2*h)/(h2+1+2*h)*(x-h)+(-1/2*cos(h)-cos(h)/(h2+1+2*h)+(sin(h)*h+sin(h)+2*cos(h)/(1+h)/(h2+1+2*h)*(2+2*h)/(h2+1+2*h)*(x-h)2便于理解，以上结果的更直观表达式可写为：【例3-15】 试对进行Taylor幂级数展开，并且观察不同阶次近似下的效果。>> x_num=-2*pi:0.1:2*pi;>> y_num=cos(x_num);>> syms x;>> y=cos(x);>> plot(x_num,y_n

23、um),axis(-2.1*pi,2.1*pi,-2,2);hold on>>i=1;for m=8:4:16 p=taylor(y,x,m);y1(i,:)=subs(p,x,x_num);i=i+1;end>>plot(x_num,y1(1,:),o);hold on>>plot(x_num,y1(2,:),s);hold on>>plot(x_num,y1(3,:),+);图 3-8 余弦函数的Taylor幂级数近似比较如函数为含有n个变量，用Taylor幂级数展开可写为：其中，为Taylor幂级数展开的中心点。Matlab中来求多变量函数

24、的Taylor幂级数展开，可以调用Maple语言中的mtaylor( )函数来求多变量函数的Taylor幂级数展开。该函数的调用格式为值得注意的是,(1)该函数调用时自变量的引号不能省略；(2)展开的最高阶次为，为原多变量函数。【例3-16】 对二元函数，求其各类Taylor幂级数展开。>> syms x y;>> fun=(y2-3*y)*exp(-(x2+y2+3*x*y);>>FUN=maple('mtaylor',fun, ' x,y ',5)FUN = -3*y+y2+3*y*x2+9*y2*x+3*y3-y2*x2

25、-3*y3*x-y4便于识别，写成一般的数学表达式为：如果求处的幂级数展开，则指令为>> syms c;>>FUN=maple('mtaylor',fun, ' x=c,y=1 ',3)结果的形式比较复杂，我们直接写为数学表达式其实，Maple中的mtaylor( )函数同样可以用于求单变量x或者y的Taylor幂级数展开。调用格式为>>FUN=maple('mtaylor',fun, ' y=c ',3)FUN =(c2-3*c)*exp(-x2-c2-3*c*x)+(c2-3*c)*exp(

26、-x2-c2-3*c*x)*(-2*c-3*x)+(2*c-3)*exp(-x2-c2-3*c*x)*(y-c)+(c2-3*c)*exp(-x2-c2-3*c*x)*(-1+2*c2+6*c*x+9/2*x2)+exp(-x2-c2-3*c*x)+(2*c-3)*exp(-x2-c2-3*c*x)*(-2*c-3*x)*(y-c)2写成数学表达式为：3.2.2 Fourier 级数展开在实际问题中，为深入研究周期函数，我们可以通过将周期函数展开成为由简单的周期函数例如三角函数组成的级数。具体的说，将周期为的周期函数用一系列以为周期的正弦函数组成的级数来表示，记为其中，都是常数。为讨论方便，将

27、正弦函数按三角公式变形，得最终可得到如果其中，并且他们都存在的话，上式就叫做的Fourier级数。通常周期函数的周期并不一定都是，那么更一般情况，周期为的函数的Fourier级数为其中，尽管Matlab和Maple都未直接提供进行Fourier系数与级数的现成函数，我们可以很容易的写出解析或数值的Fourier级数求解函数。function A,B,F=ch32fourier_series(f,x,n,lowerlimit,toplimit) if nargin=3; lowerlimit=-pi; toplimit=pi; end T=(toplimit-lowerlimit)/2; if

28、lowerlimit+toplimit f=subs(f,x,x+T+lowerlimit); end A=int(f,x,-T,T)/T;B=; F=A/2; for i=1:n An=int(f*cos(i*pi*x/T),x,-T,T)/T; Bn=int(f*sin(i*pi*x/T),x,-T,T)/T; A=A,An; B=B,Bn; F=F+An*cos(i*pi*x/T)+Bn*sin(i*pi*x/T); end if lowerlimit+toplimit F=subs(F,x,x-T-lowerlimit); end 如调用该函数，格式为其中，为给定函数，为自变量，为展开

29、项数，为的区间，如果省略，那么系统默认为。【例3-17】求函数，的Fourier级数展开。>> syms x;>> A,B,F=ch32fourier_series(x-pi/2)*(x-2*pi),x,5,0,2*pi); 薛定宇，陈阳泉。高等应用数学问题的MATLAB求解。清华大学出版社，北京，2004。这样，通过调用上面所定义的函数ch32fourier_series( )，就得到了原函数的Fourier级数展开【例3-18】考虑信号，试对该信号进行Fourier级数拟合，并观察结果。首先定义M函数function M=ch32M_fun(x,L,p)n=size

30、(x'); for i=1:n if x(i)>=0&x(i)<L/2 M(i)=p*x(i)/2;else M(i)=p*(L-x(i)/2; endend绘制出理论值>> L=2*pi;p=0.6;xx=0:2*pi/400:2*pi;>>MM=ch32M_fun(xx,L,p);>>plot(xx,MM)>>hold on>> syms x_str>> for m=1:10 n=size(xx');for i=1:nifxx(i)>=0&xx(i)<L/2a,b

31、,fun1=ch32fourier_series(p*x_str/2,x_str,m,0,L); yy(i)=subs(fun1,x_str,xx(i);elseif xx(i)>=L/2&xx(i)<=La,b,fun2=ch32fourier_series(p*(L-x_str)/2,x_str,m,0,L); yy(i)=subs(fun2,x_str,xx(i);endendplot(xx,yy);end图3-9 例3-18的结果曲线3.2.3级数求和的计算对于通式为的级数，其前项和的数学表示方法为MATLAB在符号运算工具箱中提供的求已知通项的有穷或无穷级数和的指

32、令为symsum( )，该函数的调用格式为其中，为求和级数的通项，为级数通项的自变量，和为级数求和的起始项与终止项。下面就举例说明该函数的使用。【例3-19】 求解无穷级数的和>> syms n;>> format long;>> S=double(symsum(-1)(n-1)*exp(n+1)/(2*nn+1),n,1,inf)S = 0.98906172896076这样就很容易的得到了该级数的解。【例3-20】 求解无穷级数的和>> syms n;>> S=symsum(4/(3*n-4)*(3*n+2),n,1,inf)S =

33、 -1/3再用数值方法求解该级数和，取前10000000项：>> m=1:10000000;>> S1=sum(4./(3*m-4).*(3*m+2);>> format long>> S1S1 = -0.33333337777791通过对比两种方法的结果不难看出，数值方法与解析解间存在很大的差异。即使选取极多的项数，所得到的结果依然有较大误差，达到级。尽管从表面上看累加的结果误差不会太大，事实上，由于双精度数值的有效数位数有限，所以计算通项时有效数位以后的数字加到累加量上就被忽略了，这种“大树吃小数”现象在数值分析中经常遇到，很多时候病态矩阵的

34、产生也是由于该现象，故采用纯数值方法，即使取再多项数也不能精确的得到正确的结果。以上给出的例子中的通项都是常数，直接采用累加的方式就可以近似求出结果。下面就给出一个通项中含有变量的函数级数，此问题用数值运算无法求解，必须采用符号运算工具箱来解决。【例3-21】 求解含有变量的无穷级数的和给出的求和命令如下：>> syms x n;>> f=-1/2*sin(n*x)/2n;>> H=symsum(f,n,0,inf) H =sin(x)/(4*cos(x)-5)试求下面级数与极限综合问题>> syms m n;>> limit(sym

35、sum(-12*(-1)m/m2,m,1,n),n,inf) ans = 12*hypergeom(1, 1, 1,2, 2,-1)用matlab精确显示这个表达式的值出来。 >> vpa(ans,75) ans =对比一下，直接用符号运算的无穷级数求和观察一下该问题的值。>> symsum(-12*(-1)n/n2,n,1,in)ans = pi2观察一下ans的数值结果：>> vpa(ans,75)ans =可以说这两种计算方法所得结果在数值上几乎完全一致，只是在结果的表达方式上存在一些差别。对于差别的来源，我们是这样理解的：有限项的和是无法用初等函数表

36、达的,出现了很复杂的式子,再求极限就无法得到最简表达式了。 3.3 数值微分问题运用diff( )指令求解函数的各阶导数的解析解非常方便，快捷。值得注意的是，这些函数的表达式或者说解析形式是已知的。如果函数的表达式并不事先知晓，给出的只是一些试验数据，此时要求解导数问题，运用前面的方法就无法获得问题的解析解，这样一来引入数值算法对于此类问题是非常必要的。3.3.1数值微分算法对于得到的一组时间间隔为的试验数据，当，那么某点处的导数等于相邻两点的差值除以，如果引入前向微分公式，可写为：同样地，亦可引入后向差分公式：这两种微分算法的精度都是级的。由于稍大产生的误差很大，事实证明，基于前向何后向差分

37、的数值微分算法求取高阶微分的精度较低。下面就分别来介绍具有和的中心差分算法。首先定义记为函数的形式使用Taylor级数展开，可将上式写成如下形式它的精度为，该算法的高阶微分公式为这里再给出一种具有阶精度的中心差分算法，其格式为：3.3.2中心差分方法对于中心差分方法，matlab中并没有提供现有的函数。我们可用上节中给出的具有阶精度的中心差分算法，编写出一个matlab函数，它在即使不趋于零时，依然可得到很好的近似结果。function dx,dy=chap33diffcentral(y,delt_t,n)yi_1=y 0 0 0 0 0;yi_2=0 y 0 0 0 0;yi_3=0 0 y

38、 0 0 0;yi_4=0 0 0 y 0 0;yi_5=0 0 0 0 y 0;yi_6=0 0 0 0 0 y;if n=1 dy=(-diff(yi_1)+7*diff(yi_2)+7*diff(yi_3)-diff(yi_4)/(12*delt_t); L=3; elseif n=2 dy=(-diff(yi_1)+15*diff(yi_2)-15*diff(yi_3)+diff(yi_4)/(12*delt_t2); L=3; elseif n=3 dy=(-diff(yi_1)+7*diff(yi_2)-6*diff(yi_3)-6*diff(yi_4)+7*diff(yi_5)-

39、diff(yi_6)/(8*delt_t3);L=5;elseif n=4dy=(-diff(yi_1)+11*diff(yi_2)-28*diff(yi_3)+28*diff(yi_4)-11*diff(yi_5)+diff(yi_6)/(6*delt_t4)L=5; enddy=dy(L+1:end-L);dx=(1:length(dy)+L-2-(n>2)*delt_t;【例3-22】 调用函数chap33diffcentral( )，用数值积分求函数的1-4阶导数，并且与解析解对比。对于该函数，可以很容易的得到解析解，通过生成的横坐标的向量x，通过解的结息表达式，可以很便捷的得到

40、精确解。%初始化各参数clearn=5;m=0.8;delt_x=0.05;x=0:delt_x:10;syms x1;%求解析解，同时得到各阶精确解y1=exp(x1m)*sin(n*x1);diff_y1=diff(y1);Num_diff_y1=subs(diff_y1,x1,x);diff_y2=diff(diff_y1);Num_diff_y2=subs(diff_y2,x1,x);diff_y3=diff(diff_y2);Num_diff_y3=subs(diff_y3,x1,x);diff_y4=diff(diff_y3);Num_diff_y4=subs(diff_y4,x1

41、,x);%调用本节中提供的函数求数值解，并且通过绘图比较，结果见图3-10。y=exp(x.m).*sin(n*x);dx1,y1=chap33diffcentral(y,delt_x,1);subplot(2,2,1),plot(x,Num_diff_y1,dx1,y1,':');dx2,y2=chap33diffcentral(y,delt_x,2);subplot(2,2,2),plot(x,Num_diff_y2,dx2,y2,':');dx3,y3=chap33diffcentral(y,delt_x,3);subplot(2,2,3),plot(x,

42、Num_diff_y3,dx3,y3,':');dx4,y4=chap33diffcentral(y,delt_x,4);subplot(2,2,4),plot(x,Num_diff_y4,dx4,y4,':');图 3-10 各阶导数的比较3.3.3二元函数的梯度计算给定二元函数的函数值矩阵z，则可以由matlab提供的指令gradient( )求取二元函数的梯度，z为网格数据。其调用格式为考虑坐标的情况，如果要求取的真正梯度还需要涉及生成网格的步长。假设矩阵z是建立在等间距的形式生成网格基础上的，可用下式求得的实际梯度值其中，和分别为生成网格的步长或者间距。

43、【例3-23】 考虑函数，网格数据已经给出，试由这些数据用数值方法求得梯度值，并观察误差。>> x,y=meshgrid(-3:0.1:3,-2:0.1:2);>> fun=(y.2-3*y).*exp(-(x.2+y.2+3*x.*y);>> fun_x,fun_y=gradient(fun);>> fun_x=fun_x/0.1;fun_y=fun_y/0.1;>> contour(x,y,fun,30)>> hold on;>> quiver(x,y,fun_x,fun_y)本例中给出的是带有等值线的引力

44、线图，如图3-11图 3-11 带有等值线的引力线图首先看该函数的偏导数>> syms x1 y1>>fun1=(y12-3*y1).*exp(-(x12+y12+3*x1*y1);>> fun1_x1=diff(fun1,x1)>> simple(fun1_x1)ans = -y1*(y1-3)*(2*x1+3*y1)*exp(-x12-y12-3*x1*y1)这样得到的是关于x1的导数，再看关于y1的导数。>> fun1_y1=diff(fun1,y1)>> simple(fun1_y1)ans = -exp(-x12

45、-y12-3*x1*y1)*(-2*y1+3+2*y13+3*y12*x1-6*y12-9*x1*y1)接下来要给出的是误差曲面，如图3-12和图3-13所示，程序如下，>> fun_x_exact=-y.*(y-3).*(2*x+3*y).*exp(-x.2-y.2-3*x.*y);>> fun_y_exact=-exp(-x.2-y.2-3*x.*y).*(-2*y+3+2*y.3+3*y.2.*x-6*y.2-9*x.*y);>> surf(x,y,abs(fun_x-fun_x_exact)>> figure;>> surf(

46、x,y,abs(fun_y-fun_y_exact)图 3-12误差曲面图 3-13的误差曲面下面我们再观察一下如果把网格加密，误差会如何变化。网格加密后的误差曲面，如图3-14和图3-15所示。>> x,y=meshgrid(-3:0.05:3,-2:0.05:2);>> fun=(y.2-3*y).*exp(-(x.2+y.2+3*x.*y);>> fun_x,fun_y=gradient(fun);>> fun_x=fun_x/0.05;fun_y=fun_y/0.05;>> fun_x_exact=-y.*(y-3).*(2*

47、x+3*y).*exp(-x.2-y.2-3*x.*y);>> fun_y_exact=-exp(-x.2-y.2-3*x.*y).*(-2*y+3+2*y.3+3*y.2.*x-6*y.2-9*x.*y);>> figure>> surf(x,y,abs(fun_x-fun_x_exact); axis(-4 4 -2 2 0 150)>> figure;>> surf(x,y,abs(fun_y-fun_y_exact); axis(-4 4 -2 2 0 2000)图 3-14 误差曲面图 图 3-15 误差曲面可以观察到，通过

48、加密网格，误差比未加密前减小了很多，通过最大误差的变化就可以看出来。3.4 数值积分问题.数值积分主要用于计算无法解析求解的丁积分的近似解，是工程师和科学家使用的基本工具。例如，在概率论中的正态分布，它的概率密度函数为，其中为常数，其概率分布函数为由于不存在解析表达式，一般通过数值方法来得到其近似值。如，在区间，时，的数值近似为：观察一下在区间上围成的图线，如图3-16。图 3-16 在区间上围成的图线表3-1列出了在区间0，5的几个点上通过数值积分得到的一些近似值。表3-1 数值积分的近似值x3.003.50.226627352376868213484.00.3085375387259869

49、15634.505.00.500000000000000031235.50.598706325682923761646.006.50.773372647623131848987.007.50.894350226333144798158.003.4.1由给定数据进行梯形求积积分的目的是，通过在有限个采样点上计算的值来逼近在区间上的定积分，一元函数定积分的数学表示为如果被积函数理论上不可积，无法得到该积分的解析解，此时就要求助于数值方法来求解。求解定积分的数值方法有很多种，包括梯形公式、辛普森公式、龙贝格积分、自适应积分和高斯-勒让德积分等。这些方法都有一个共同的基本思路，都是将整个积分空间分割成

50、若干个子空间，其中，这样一来整个积分问题就可以归结为下面的求和形式。经过这样的变换，每个小的子空间上的值都可以近似地求解出来，采用梯形近似的方法求每一个小子空间的积分方法是最方便简捷的。Matlab中实现梯形积分算法很容易。设有向量，可用以下语句得到积分值另外，matlab中提供的指令trapz( )同样可直接用梯形法求解积分问题，该函数的调用格式为其中，为维数相同的行向量或列向量。下面就通过算例来简单演示一下运用以上两个命令求积分。【例3-24】 用梯形法求函数在区间区间内，求函数的定积分值。>> x_num=1:0.01:2'>>y=exp(x_num) x

51、_num.2 sin(x_num);>>x=x_num x_num x_num;>>I=sum(2*y(1:end-1,:)+diff(y).*diff(x)/2I =4.67081319352565 2.33335000000000 0.95644117199248用指令trapz( )作对比>>I1=trapz(x_num,y)I1 = 4.67081319352565 2.33335000000000 0.95644117199248可注意到，结果完全一致。【例3-25】 用定步长法求解积分首先观察一下该函数的曲线，见图3-17，>>x=0

52、:pi/2/100:pi/2;>> y=exp(1.5*x).*sin(30*x);>> plot(x,y)图 3-17 被积函数的曲线在不同步长，采用下面的语句得到近似积分结果，并且通过与精确解的对比观察它们的积分精度，>> syms x;>> S=int(exp(3/2*x)*sin(30*x),0,pi/2); S = 40/1203*exp(pi)(3/4)+40/1203S就是积分的精确解>>S_err=;>>delt_x=pi/2*0.1 0.01 0.001 0.0001 0.00001 0.000001;&

53、gt;>for i=1:6 x1=0:delt_x(i):pi/2; y=exp(3/2*x1).*sin(30*x1); S_app(i)=trapz(x1,y);end>>S_err=S-S_app;>>vpa(S_err,10)ans = 1.266643304,.71513767e-2, .712531e-4, .7124e-6, .70e-8, -.1e-9很明显，随着步长的减小，计算精度逐渐增加，通过误差的变化就可以看出来。当然，随着步长的减小，计算时间也在逐渐的增加，如果想进一步提高计算精度，势必要增加计算时间。3.4.2单变量数值积分求解对于单变量

54、函数的数值积分问题可以采用一般数值积分中的算法进行求解，如采用辛普森公式进行逼近积分。一般的用下面公式逼近 其中。 在MATLAB中来实现该函数function Z=chap34simpson_rule(f,a,b,tol0)h=(b-a)/2;C=zeros(1,3);C=feval(f,a (a+b)/2 b);S=h*(C(1)+4*C(2)+C(3)/3;S2=S;tol1=tol0;err=tol0;Z=a b S S2 errtol1;下面举例来说明该函数chap34simpson_rule( )的应用。【例3-26】考虑函数，其解析不可积，需采用数值方法求解。对于如何描述被积函数

55、，通常有两种方法，第一建立一个matlab函数并将其存为文件，第二就是用inline()指令来定义被积函数。如果采用前者，那么将下列内容存成文件function y=chap34erf_define(x)y=2/sqrt(pi)*exp(-x.2);这样将含有上述内容的文件存入一个chap34erf_define.m文件。如果用第二种方法，可用下面的语句来定义被积函数。>>f=inline('2/sqrt(pi)*exp(-x.2)','x');这样无需建立一个单独的matlab文件，比前者更加便捷。接下来就来求解该问题>>f=inlin

56、e('2/sqrt(pi)*exp(-x.2)','x');>>Z=chap34simpson_rule(f,0,1.5,1e-6);>> Z(3)ans =9.547584332749929e-001为比较结果，用MATLAB指令quad( )来观察结果，其调用格式为其中，为被积函数，为上下限，为限定精度。>> y = quad(f, 0, 1.5, 1e-10)y =可以观察到误差为1.174%，通常情况下可以接受。还可以用以下公式来进行逼近积分这样根据上面公式就可以给出辛普森公式的自适应积分找参考文献function SRmat,quad,err=chap34simpson_adapt(f,a,b,tol)%初始化SRmat = zeros(30,6);iterating=0;done=1;SRvec