1、微积分的计算,微积分的符号运算,Matlab定义符号运算和数值运算。如果运算过程中有符号变量，则运算过程为符号运算。符号变量可以利用两种方法定义。,(1) 利用函数sym定义符号表达式语句格式 变量=sym(字符串)sym函数可以定义符号变量、符号常数和一般的符号表达式。, x=sym(1/3) x =1/3 y=x+1/2 %符号运算 y =5/6,x=1/3x= 0.3333 y=1/2+1/3 %数值计算y = 0.8333,例1：比较符号运算和数值计算,例2：分别利用符号运算和数值运算计算如下的和 并比较其计算速度。,注：Matlab语言中利用tic 和toc 计算时间，语言格式为 t
2、ic 程序块 toc显示器显示tic和toc之间的程序块运行时间,符号运算程序：s=sym(0);ticfor k=1:1000 s=s+1/k;endtocs,运行结果：elapsed time is 17.471170 seconds. s = 5336291328229478504559104562404298040965247228038426009710134924845626888949710175750609790198503569140908873155046809837844217211788500946430234432656602250210027842563285208
3、1405544941210442510142672770294774712708917963967779610453224692426866468888281582071984897105110796873249319155529397017508931564519976085734473014183284011724412280649074307703736683170055800293659235088589360235285852808160759574737836655413175508131522517/7128865274665093053166384155714272920668
4、358861885893040452001991154324087581111499476444151913871586911717817019575256512980264067621009251465871004305131072686268143200196609974862745937188343705015434452523739745298963145674982128236956232823794011068809262317708861979540791247754558049326475737829923352751796735248042463638051137034331
5、214781746850878453485678021888075373249921995672056932029099390891687487672697950931603520000,数值运算程序 s=0; tic for k=1:1000 s=s+1/k; end toc s运行结果：Elapsed time is 0.000015 seconds.s =7.4855,可以看出，数值计算的计算速度远远高于符号运算。另外，符号运算只能计算非常简单的问题，复杂问题只能利用数值方法。但由于符号运算得到的是解析式，在许多情况下有利于问题的进一步分析。,利用语句 syms 定义符号变量语句形式 s
6、yms x y zsyms语句 一次可定义多个符号变量，但不能定义符号常数和表达式,微积分计算问题的符号运算,极限运算语句格式 G= limit(F); %表达式F中变量趋向于0的极限 G= limit(F,a); %表达式F中变量趋向于a的极限 G= limit(F,v,a); %表达式F中变量v趋向于a的极限 G= limit(F,v,a,right); %表达式F中变量趋向于a的右极限。,例3：计算 syms xf=(1-x)/(1+x)(2/x) f =(-(x - 1)/(x + 1)(2/x) limit(f,0) ans =1/exp(4),例4：计算 syms x y z f=
7、sin(x*y)/sin(x*z) f =sin(x*y)/sin(x*z) limit(f,x,0) ans =y/z,导数求导数的语句为 df=diff(F) %求表达式F的一阶导数 df=diff(F,n) %求表达式F的n阶导数 df=diff(F,u,n) %求表达式F的关于变量u的n阶导数,例5：求函数 f=sin(x)/x 的二阶导数。 f=sin(x)/xf = sin(x)/x diff(f,2)ans =(2*sin(x)/x3 - sin(x)/x - (2*cos(x)/x2,例6:求函数 关于x和y的二阶偏导数。f=(x-y)/(x+y)f = (x - y)/(x
8、+ y) dfxy=diff(diff(f,y),x)dfxy =(2*(x - y)/(x + y)3,注：如果表达式F有多个符号变量，且没有指明对哪一个变量求导，则Matlab规定对离x最近的变量求导。其他运算也有类似规定。,例7：求 f=x2*sin(x) exp(x)*cos(x) log(1+x2) atan(x); diff(f)ans = x2*cos(x) + 2*x*sin(x), exp(x)*cos(x) - exp(x)*sin(x) (2*x)/(x2 + 1), a/(x2 + 1),Matlab的微积分符号运算都可以对数组进行。,函数的积分积分符号运算的基本语句
9、int(F)； %求函数表达式F的不定积分 int(F,v)； %求函数表达式F关于变量v的不定积分 int(F,a,b)； %求函数表达式F在区间a,b上的定积分 int(F,v,a,b); %求函数表达式F在区间a,b上的关于变量v的定积分,例8：求不定积分 syms a b x I=int(exp(a*x)*cos(b*x),x)I = (exp(a*x)*(a*cos(b*x) + b*sin(b*x)/(a2 + b2),例9：求定积分 syms x I=int(sin(x)15,x,0,pi/2) I =2048/6435,对含参数的运算，符号计算的优势在于可以提供一个用于进一步分
10、析的解析表达式。,级数的符号计算级数符号计算的基本语句 S=symsum(g)； % 求通项g的无穷和 S=symsum(g,n)； % 求通项g的对整变量n的无穷和 S=symsum(g,a,b)； % 求通项g的下标在a和b之间的项的和 S=symsum(g,n,a,b)； % 求通项g的下标变量n在a和b之间的项的和,例10:求级数的和 syms n symsum(1/n2,1,inf) ans = pi2/6,例11: 求幂级数的和 syms x n f=xn/n f = xn/n S=symsum(f,n,1,inf) S = piecewise(1 1, -log(1 - x),T
11、aylor展式Taylor展式的基本语句 F=taylor(f); %求函数f 的马克劳林级数的前6项 F=taylor(f,n); %求函数f 的马克劳林级数的前n项 F=taylor(f,n,x0); %求函数f 的在x0展开的Taylor级数的前n项 F=taylor(f,v,n,v0); %求函数f对变量v在v0展开的Taylor级数的前n项,例12: 求函数 的在0点的taylor展式前8项 syms x S=taylor(1-x)/(1+x),8) S = 2*x6 - 2*x7 - 2*x5 + 2*x4 - 2*x3 + 2*x2 - 2*x + 1,微积分的数值计算,符号运算
12、有很大的局限性。第一是计算是纯推理过程，效率非常低。第二是大量的问题无法利用解析的形式给出，例如我们经常使用的概率积分 就无法利用牛顿-莱布尼兹公式推导得到。,相比之下，数值计算在实际计算中有符号运算无法比拟的优点。一是效率高，对公式稍微复杂一些或规模稍微大一些的问题，只能利用数值计算。二是我们经常遇到的计算问题，都可以利用数值计算得到满意的解，而符号运算只能处理其中一小部分问题。,数值计算需要对连续变量离散求解！数值计算需要设计数值算法,数值导数,数值计算是在离散节点上的近似计算，通常对连续的区间要首先进行离散。对计算区间a.b上的导数。,首先将区间离散成等距的节点其中步长 与下标k无关。,
13、一阶和二阶导数的计算由公式 给出, 误差为 O(h2)。,Matlab的数值导数计算可以利用语句diff实现。diff语句对非符号计算表示向前差分，即对 函数diff(x)输出diff(x,2)输出 diff(diff(x),即 一般, diff(x,k) 的计算可依次类推。,注： diff(x) 的元素个数比x少 1！,一阶导数的计算一阶导数可以利用 计算。由公式df的第i个元素与 f(x) 在 xi+1/2点的值最接近，误差是O(h2),例13：利用diff语句计算函数 y=sin(x) 的导数，并计算误差。参考程序：x=0:pi/100:2*pi;y=sin(x);x1=pi/200:p
14、i/100:2*pi;y1=diff(y)./diff(x);y11=sin(x1);subplot(1,2,1)plot(x1,y1,x1,y11)subplot(1,2,2)y12=cos(x1);plot(x1,y1-y12)legend(df,f)legend(error) 结果图：,二阶导数的计算语句 输出二阶向前差分，即 可以利用 来计算y(x)的二阶导数。,计算的是n-2个内点上的二阶导数，不包含两个端点上的导数！,例14：计算 y=sin(x)的二阶导数，并与 sin(x)比较。参考程序：x=0:pi/100:2*pi;y=sin(x);dx=diff(x);d2y=diff(
15、y,2)./dx(2:end).2;subplot(1,2,1)plot(x,y,x(2:end-1),d2y)legend(func,numerd2f)subplot(1,2,2)plot(x(2:end-1),d2y+sin(x(2:end-1)legend(error),数值积分一元数值积分Matlab中的数值积分函数quad采用自适应复化Simpson公式，语句格式为 I=quad(fun,a,b) %计算函数fun在a,b上的定积分。 I=quad(fun,a,b,tol) %计算函数fun在a,b上精度tol的定积分。,fun的 取法： 过去的取法是一个字符串作为所调用的函数名字。
16、如 quad(sin, 0, pi/2) 在matlab7以上的版本中引入了函数句柄，可以在函数quad中利用指示一个函数句柄。例如 sin %调用以sin为函数名的函数子程序 (x,y)(x.3+y.3) %建立并调用一个函数句柄 (x)myfun(x,a) %调用含有参数的子程序myfun,例15：计算积分 k=1,2,3,4,5并与精确值比较。,I=zeros(1,5);for k=1:5 g=strcat(sin(x).,num2str(k); I(k)=quad(g,0,pi/2);endI0=1,pi/4,2/3,3*pi/16,8/15;subplot(1,2,1)plot(I-
17、I0,p)legend(errors of quad) for k=1:5 g=strcat(sin(x).,num2str(k); I(k)=quad(g,0,pi/2,1e-12); end subplot(1,2,2)plot(I-I0,p)legend(errors of quad),积分函数quadv积分函数quadv类似于函数quad，但它可以对含有向量参数的函数积分，而且可以对矩阵函数积分。,例16：计算积分k=1,2,3,4,5,这里k 可以看做参数向量!,参考语句： format rat I=quadv(sin(x).1:5,0,pi/2)I = 1 355/452 2/3
18、1065/1808 8/15,例17：计算矩阵函数 的积分,参考语句： format rat I=quadv(x)exp(-x).*cos(x),(1-x)./sqrt(1+x.2);x.*log(1+x),sin(x)./x,1,2)I = 249/7261 -548/2111 1075/769 689/1045,可用于广义积分的积分函数quadgk函数quadgk采用Gauss-Kronrod数值积分方法，可以解决无穷区间、无穷间断点等问题。其使用方法与quad类似，但积分限可以取inf,例18：计算 参考语句： I=quadgk(x)1./(3*x.2+2),1,inf)I =0.279
19、5,多重积分Matlab的多重积分函数为 I=dblquad(fun,xmin,xmax,ymin,ymax) %矩形区域的二重积分 I=dblquad(fun,xmin,xmax,ymin,ymax,tol) %矩形区域的二重积分，误差限tol I=triplequad(fun, xmin,xmax,ymin,ymax,zmin,zmax) %矩形区域的三重积分I=triplequad(fun, xmin,xmax,ymin,ymax,zmin,zmax),例19：计算下述二重积分 （1） （2）, I1=dblquad(x,y)(x.*exp(y)-y.*exp(x),0,1,1,3)I1 = 1.8105 I2=dblquad(x,y)(2*x.2+y)./sqrt(1+x.2+y.2).*(x.2+y.2 Q1=triplequad(x,y,z)(1-x-2*y+3*z),0,1,0,1,0,1)Q1 = 1 Q2=triplequad(x,y,z)(x+y+z)./sqrt(1+x.2+y.2+z.2).*(x.2+y.2+z.2=1),-1,1,-1,1,0,1)Q2 =0.6134,
