數(shù)學(xué)實(shí)驗(yàn)微積分問(wèn)題的計(jì)算機(jī)求解PPT課件

1、3.1 3.1 微積分問(wèn)題的解析解 極限問(wèn)題的解析解 單變量函數(shù)的極限 格式1： L= limit( fun, x, x0) 格式2： L= limit( fun, x, x0, left 或 right)第1頁(yè)/共103頁(yè) 例: 試求解極限問(wèn)題 syms x a b; f=x*(1+a/x)x*sin(b/x); L=limit(f,x,inf) L = exp(a)*b 例:求解單邊極限問(wèn)題 syms x; limit(exp(x3)-1)/(1-cos(sqrt(x-sin(x),x,0,right) ans =12第2頁(yè)/共103頁(yè) 在(-0.1,0.1)區(qū)間繪制出函數(shù)曲線： x=-0

2、.1:0.001:0.1; y=(exp(x.3)-1)./(1-cos(sqrt(x-sin(x);Warning: Divide by zero.(Type warning off MATLAB:divideByZero to suppress this warning.) plot(x,y,-,0,12,o)第3頁(yè)/共103頁(yè)多變量函數(shù)的極限： 格式： L1=limit(limit(f,x,x0),y,y0) 或 L1=limit(limit(f,y,y0), x,x0) 如果x0 或y0不是確定的值，而是另一個(gè)變量的函數(shù)，如x-g(y)，則上述的極限求取順序不能交換。第4頁(yè)/共103頁(yè)

3、 例：求出二元函數(shù)極限值 syms x y a; f=exp(-1/(y2+x2) *sin(x)2/x2*(1+1/y2)(x+a2*y2); L=limit(limit(f,x,1/sqrt(y),y,inf)L =exp(a2)第5頁(yè)/共103頁(yè)函數(shù)導(dǎo)數(shù)的解析解 函數(shù)的導(dǎo)數(shù)和高階導(dǎo)數(shù) 格式： y=diff(fun,x) %求導(dǎo)數(shù)(默認(rèn)為1階) y= diff(fun,x,n) %求n階導(dǎo)數(shù) 例： 一階導(dǎo)數(shù)： syms x; f=sin(x)/(x2+4*x+3); f1=diff(f); pretty(f1)第6頁(yè)/共103頁(yè) cos(x) sin(x) (2 x + 4) - - -

4、2 2 2 x + 4 x + 3 (x + 4 x + 3)原函數(shù)及一階導(dǎo)數(shù)圖： x1=0:.01:5; y=subs(f, x, x1); y1=subs(f1, x, x1); plot(x1,y,x1,y1,:)更高階導(dǎo)數(shù)： tic, diff(f,x,100); tocelapsed_time = 4.6860第7頁(yè)/共103頁(yè) 原函數(shù)4階導(dǎo)數(shù) f4=diff(f,x,4); pretty(f4) 2 sin(x) cos(x) (2 x + 4) sin(x) (2 x + 4) - + 4 - - 12 - 2 2 2 2 3 x + 4 x + 3 (x + 4 x + 3)

5、(x + 4 x + 3) 3 sin(x) cos(x) (2 x + 4) cos(x) (2 x + 4) + 12 - - 24 - + 48 - 2 2 2 4 2 3 (x + 4 x + 3) (x + 4 x + 3) (x + 4 x + 3) 4 2 sin(x) (2 x + 4) sin(x) (2 x + 4) sin(x) + 24 - - 72 - + 24 - 2 5 2 4 2 3 (x + 4 x + 3) (x + 4 x + 3) (x + 4 x + 3)第8頁(yè)/共103頁(yè) 多元函數(shù)的偏導(dǎo)： 格式： f=diff(diff(f,x,m),y,n) 或

6、f=diff(diff(f,y,n),x,m) 例： 求其偏導(dǎo)數(shù)并用圖表示。 syms x y z=(x2-2*x)*exp(-x2-y2-x*y); zx=simple(diff(z,x)zx = -exp(-x2-y2-x*y)*(-2*x+2+2*x3+x2*y-4*x2-2*x*y)第9頁(yè)/共103頁(yè) zy=diff(z,y)zy =(x2-2*x)*(-2*y-x)*exp(-x2-y2-x*y) 直接繪制三維曲面 x,y=meshgrid(-3:.2:3,-2:.2:2); z=(x.2-2*x).*exp(-x.2-y.2-x.*y); surf(x,y,z), axis(-3

7、3 -2 2 -0.7 1.5) 第10頁(yè)/共103頁(yè) contour(x,y,z,30), hold on % 繪制等值線 zx=-exp(-x.2-y.2-x.*y).*(-2*x+2+2*x.3+x.2.*y-4*x.2-2*x.*y); zy=-x.*(x-2).*(2*y+x).*exp(-x.2-y.2-x.*y); % 偏導(dǎo)的數(shù)值解 quiver(x,y,zx,zy) % 繪制引力線第11頁(yè)/共103頁(yè) 例 syms x y z; f=sin(x2*y)*exp(-x2*y-z2); df=diff(diff(diff(f,x,2),y),z); df=simple(df); p

8、retty(df) 2 2 2 2 2 -4 z exp(-x y - z ) (cos(x y) - 10 cos(x y) y x + 4 2 4 2 2 4 2 2sin(x y) x y+ 4 cos(x y) x y - sin(x y)第12頁(yè)/共103頁(yè) 多元函數(shù)的Jacobi矩陣: 格式：J=jacobian(Y,X)其中，X是自變量構(gòu)成的向量，Y是由各個(gè)函數(shù)構(gòu)成的向量。第13頁(yè)/共103頁(yè) 例：試推導(dǎo)其 Jacobi 矩陣 syms r theta phi; x=r*sin(theta)*cos(phi); y=r*sin(theta)*sin(phi); z=r*cos(t

9、heta); J=jacobian(x; y; z,r theta phi) J = sin(theta)*cos(phi), r*cos(theta)*cos(phi), -r*sin(theta)*sin(phi) sin(theta)*sin(phi), r*cos(theta)*sin(phi), r*sin(theta)*cos(phi) cos(theta), -r*sin(theta), 0 第14頁(yè)/共103頁(yè) 隱函數(shù)的偏導(dǎo)數(shù)： 格式：F=-diff(f,xj)/diff(f,xi)第15頁(yè)/共103頁(yè) 例： syms x y; f=(x2-2*x)*exp(-x2-y2-x*

10、y); pretty(-simple(diff(f,x)/diff(f,y) 3 2 2 -2 x + 2 + 2 x + x y - 4 x - 2 x y - - x (x - 2) (2 y + x)第16頁(yè)/共103頁(yè)積分問(wèn)題的解析解 不定積分的推導(dǎo)： 格式： F=int(fun,x) 例：用diff() 函數(shù)求其一階導(dǎo)數(shù)，再積分，檢驗(yàn)是否可以得出一致的結(jié)果。 syms x; y=sin(x)/(x2+4*x+3); y1=diff(y); y0=int(y1); pretty(y0) % 對(duì)導(dǎo)數(shù)積分 sin(x) sin(x) - 1/2 - + 1/2 - x + 3 x + 1第

11、17頁(yè)/共103頁(yè) 對(duì)原函數(shù)求4 4 階導(dǎo)數(shù)，再對(duì)結(jié)果進(jìn)行4 4次積分 y4=diff(y,4); y0=int(int(int(int(y4); pretty(simple(y0) sin(x) - 2 x + 4 x + 3第18頁(yè)/共103頁(yè) 例:證明 syms a x; f=simple(int(x3*cos(a*x)2,x)f = 1/16*(4*a3*x3*sin(2*a*x)+2*a4 *x4+6*a2*x2*cos(2*a*x)-6*a*x*sin(2*a*x)-3*cos(2*a*x)-3)/a4 f1=x4/8+(x3/(4*a)-3*x/(8*a3)*sin(2*a*x)

12、+. (3*x2/(8*a2)-3/(16*a4)*cos(2*a*x); simple(f-f1) % 求兩個(gè)結(jié)果的差ans = -3/16/a4第19頁(yè)/共103頁(yè) 定積分與無(wú)窮積分計(jì)算: 格式： I=int(f,x,a,b) 格式： I=int(f,x,a,inf)第20頁(yè)/共103頁(yè) 例： syms x; I1=int(exp(-x2/2),x,0,1.5) 無(wú)解I1 =1/2*erf(3/4*2(1/2)*2(1/2)*pi(1/2) vpa(I1,70) ans = I2=int(exp(-x2/2),x,0,inf) I2 =1/2*2(1/2)*pi(1/2) 2/ 2( )x

13、fxe202( )xterf xe dt第21頁(yè)/共103頁(yè) 多重積分問(wèn)題的MATLAB求解 例： syms x y z; f0=-4syms x y z; f0=-4* *z z* *exp(-x2exp(-x2* *y-y-z2)z2)* *(cos(x2(cos(x2* *y)-10y)-10* *cos(x2cos(x2* *y)y)* *y y* *x2+.x2+. 4 4* *sin(x2sin(x2* *y)y)* *x4x4* *y2+4y2+4* *cos(x2cos(x2* *y)y)* *x4x4* *y2-y2-sin(x2sin(x2* *y);y); f1=int(

14、f0,z);f1=int(f1,y);f1=int(f1,x); f1=int(f0,z);f1=int(f1,y);f1=int(f1,x); f1=simple(int(f1,x) f1=simple(int(f1,x)f1 =f1 = exp(-x2 exp(-x2* *y-z2)y-z2)* *sin(x2sin(x2* *y)y)第22頁(yè)/共103頁(yè) f2=int(f0,z); f2=int(f2,x); f2=int(f2,x); f2=simple(int(f2,y)f2 =2*exp(-x2*y-z2)*tan(1/2*x2*y)/(1+tan(1/2*x2*y)2) simp

15、le(f1-f2)ans =0 順序的改變使化簡(jiǎn)結(jié)果不同于原函數(shù)，但其誤差為0，表明二者實(shí)際完全一致。這是由于積分順序不同，得不出實(shí)際的最簡(jiǎn)形式。第23頁(yè)/共103頁(yè) 例： syms x y z int(int(int(4*x*z*exp(-x2*y-z2),x,0,1),y,0,pi),z,0,pi)ans =(Ei(1,4*pi)+log(pi)+eulergamma+2*log(2)*pi2*hypergeom(1,2,-pi2)Ei(n,z)為指數(shù)積分，無(wú)解析解，但可求其數(shù)值解： vpa(ans,60) ans = 第24頁(yè)/共103頁(yè)3.2 3.2 函數(shù)的級(jí)數(shù)展開(kāi)與 級(jí)數(shù)求和問(wèn)題求解

16、冪級(jí)數(shù)展開(kāi)級(jí)數(shù)展開(kāi)級(jí)數(shù)求和的計(jì)算第25頁(yè)/共103頁(yè)冪級(jí)數(shù)展開(kāi) 單變量函數(shù)的 TaylorTaylor 冪級(jí)數(shù)展開(kāi)第26頁(yè)/共103頁(yè)第27頁(yè)/共103頁(yè)例: syms x; f=sin(x)/(x2+4*x+3); y1=taylor(f,x,9); pretty(y1) 2 23 3 34 4 4087 5 3067 6 515273 7 386459 8 1/3 x - 4/9 x + - x - - x + -x - - x +- x - - x 54 81 9720 7290 1224720 918540第28頁(yè)/共103頁(yè) taylor(f,x,9,2)ans = syms a;

17、taylor(f,x,5,a) % 結(jié)果較冗長(zhǎng)，顯示從略ans =sin(a)/(a2+3+4*a) +(cos(a)-sin(a)/(a2+3+4*a)*(4+2*a)/(a2+3+4*a)*(x-a) +(-sin(a)/(a2+3+4*a)-1/2*sin(a)-(cos(a)*a2+3*cos(a)+4*cos(a)*a-4*sin(a)-2*sin(a)*a)/(a2+3+4*a)2*(4+2*a)/(a2+3+4*a)*(x-a)2+第29頁(yè)/共103頁(yè)例:對(duì)y=sinx進(jìn)行Taylor冪級(jí)數(shù)展開(kāi),并觀察不同階次的近似效果。 x0=-2*pi:0.01:2*pi; y0=sin(x

18、0); syms x; y=sin(x); plot(x0,y0,r-.), axis(-2*pi,2*pi,-1.5,1.5); hold on for n=8:2:16 p=taylor(y,x,n), y1=subs(p,x,x0); line(x0,y1) endp =x-1/6*x3+1/120*x5-1/5040*x7p =x-1/6*x3+1/120*x5-1/5040*x7+1/362880*x9p =x-1/6*x3+1/120*x5-1/5040*x7+1/362880*x9-1/39916800*x11p =x-1/6*x3+1/120*x5-1/5040*x7+1/36

19、2880*x9-1/39916800*x11+1/6227020800*x13 第30頁(yè)/共103頁(yè)p =第31頁(yè)/共103頁(yè)多變量函數(shù)的TaylorTaylor 冪級(jí)數(shù)展開(kāi) 多變量函數(shù) 在的Taylor冪級(jí)數(shù)的展開(kāi)12(,)nf x xx12(,)na aa第32頁(yè)/共103頁(yè) 例：？ syms x y; f=(x2-2*x)*exp(-x2-y2-x*y); F=maple(mtaylor,f,x,y,8)F = mtaylor(x2-2*x)*exp(-x2-y2-x*y),x, y,8)第33頁(yè)/共103頁(yè) maple(readlib(mtaylor);讀庫(kù)，把函數(shù)調(diào)入內(nèi)存 F=map

20、le(mtaylor,f,x,y,8) F =-2*x+x2+2*x3-x4-x5+1/2*x6+1/3*x7+2*y*x2+2*y2*x-y*x3-y2*x2-2*y*x4-3*y2*x3-2*y3*x2-y4*x+y*x5+3/2*y2*x4+y3*x3+1/2*y4*x2+y*x6+2*y2*x5+7/3*y3*x4+2*y4*x3+y5*x2+1/3*y6*x syms a; F=maple(mtaylor,f,x=1,y=a,3); F=maple(mtaylor,f,x=a,3)F =(a2-2*a)*exp(-a2-y2-a*y)+(a2-2*a)*exp(-a2-y2-a*y)

21、*(-2*a-y)+(2*a-2)*exp(-a2-y2-a*y)*(x-a)+(a2-2*a)*exp(-a2-y2-a*y)*(-1+2*a2+2*a*y+1/2*y2)+exp(-a2-y2-a*y)+(2*a-2)*exp(-a2-y2-a*y)*(-2*a-y)*(x-a)2第34頁(yè)/共103頁(yè)級(jí)數(shù)展開(kāi)第35頁(yè)/共103頁(yè)function A,B,F=fseries(f,x,n,a,b)if nargin=3, a=-pi; b=pi; endL=(b-a)/2; if a+b, f=subs(f,x,x+L+a); end變量區(qū)域互換A=int(f,x,-L,L)/L; B=; F

22、=A/2; for i=1:n an=int(f*cos(i*pi*x/L),x,-L,L)/L; bn=int(f*sin(i*pi*x/L),x,-L,L)/L; A=A, an; B=B,bn; F=F+an*cos(i*pi*x/L)+bn*sin(i*pi*x/L);endif a+b, F=subs(F,x,x-L-a); end 換回變量區(qū)域第36頁(yè)/共103頁(yè)例: syms x; f=x*(x-pi)*(x-2*pi); A,B,F=fseries(f,x,6,0,2*pi)A = 0, 0, 0, 0, 0, 0, 0 B = -12, 3/2, -4/9, 3/16, -1

23、2/125, 1/18 F =12*sin(x)+3/2*sin(2*x)+4/9*sin(3*x)+3/16*sin(4*x)+12/125*sin(5*x)+1/18*sin(6*x)第37頁(yè)/共103頁(yè)例: syms x; f=abs(x)/x; % 定義方波信號(hào) xx=-pi:pi/200:pi; xx=xx(xx=0); xx=sort(xx,-eps,eps); % 剔除零點(diǎn) yy=subs(f,x,xx); plot(xx,yy,r-.), hold on % 繪制出理論值并保持坐標(biāo)系 for n=2:20 a,b,f1=fseries(f,x,n), y1=subs(f1,x,

24、xx); plot(xx,y1)end第38頁(yè)/共103頁(yè)a = 0, 0, 0b = 4/pi, 0f1 =4/pi*sin(x)a = 0, 0, 0, 0 b = 4/pi, 0, 4/3/pif1 =4/pi*sin(x)+4/3/pi*sin(3*x)第39頁(yè)/共103頁(yè)級(jí)數(shù)求和的計(jì)算 是在符號(hào)工具箱中提供的第40頁(yè)/共103頁(yè)例:計(jì)算 format long; sum(2.0:63) %數(shù)值計(jì)算ans = 1.844674407370955e+019 sum(sym(2).0:200) % 或 syms k; symsum(2k,0,200)把2定義為符號(hào)量可使計(jì)算更精確ans =

25、 syms k; symsum(2k,0,200)ans =第41頁(yè)/共103頁(yè)例:試求解無(wú)窮級(jí)數(shù)的和 syms n; s=symsum(1/(3*n-2)*(3*n+1),n,1,inf)%采用符號(hào)運(yùn)算工具箱s =1/3 m=1:10000000; s1=sum(1./(3*m-2).*(3*m+1);%數(shù)值計(jì)算方法,雙精度有效位16,“大數(shù)吃小數(shù)”,無(wú)法精確 format long; s1 % 以長(zhǎng)型方式顯示得出的結(jié)果s1 = 0.33333332222165第42頁(yè)/共103頁(yè)例:求解 syms n x s1=symsum(2/(2*n+1)*(2*x+1)(2*n+1),n, 0,in

26、f); simple(s1) % 對(duì)結(jié)果進(jìn)行化簡(jiǎn)，MATLAB 6.5 及以前版本因本身 bug 化簡(jiǎn)很麻煩ans =log(2*x+1)2)(1/2)+1)/(2*x+1)2)(1/2)-1)%實(shí)際應(yīng)為log(x+1)/x)第43頁(yè)/共103頁(yè)例:求 syms m n; limit(symsum(1/m,m,1,n)-log(n),n,inf)ans =eulergamma vpa(ans, 70) % 顯示 70 位有效數(shù)字ans = 第44頁(yè)/共103頁(yè)3.3 3.3 數(shù)值微分 數(shù)值微分算法 向前差商公式： 向后差商公式第45頁(yè)/共103頁(yè)兩種中心公式：2342342( )( )( )/

27、 2!( )/3!()( )2( )( )( )/ 2!( )/3!()2( )( )3!f xtfxt fxt fotf xtf xtfxt fxt fotttfxf 第46頁(yè)/共103頁(yè)第47頁(yè)/共103頁(yè)第48頁(yè)/共103頁(yè)中心差分方法及其 MATLAB MATLAB 實(shí)現(xiàn) function dy,dx=diff_ctr(y, Dt, n) yx1=y 0 0 0 0 0; yx2=0 y 0 0 0 0; yx3=0 0 y 0 0 0; yx4=0 0 0 y 0 0; yx5=0 0 0 0 y 0; yx6=0 0 0 0 0 y; switch n case 1 dy = (-

28、diff(yx1)+7*diff(yx2)+7*diff(yx3)- diff(yx4)/(12*Dt); L0=3; case 2 dy=(-diff(yx1)+15*diff(yx2)- 15*diff(yx3) +diff(yx4)/(12*Dt2);L0=3; 數(shù)值計(jì)算diff(X)表示數(shù)組X相鄰兩數(shù)的差第49頁(yè)/共103頁(yè) case 3 dy=(-diff(yx1)+7*diff(yx2)-6*diff(yx3)-6*diff(yx4)+. 7*diff(yx5)-diff(yx6)/(8*Dt3); L0=5; case 4 dy = (-diff(yx1)+11*diff(yx2

29、)-28*diff(yx3)+28* diff(yx4)-11*diff(yx5)+diff(yx6)/(6*Dt4); L0=5; end dy=dy(L0+1:end-L0); dx=(1:length(dy)+L0-2-(n2)*Dt;調(diào)用格式： y為 等距實(shí)測(cè)數(shù)據(jù)， dy為得出的導(dǎo)數(shù)向量， dx為相應(yīng)的自變量向量，dy、dx的數(shù)據(jù)比y短 。,_( , )yxdddiffctr yt n第50頁(yè)/共103頁(yè) 例：求導(dǎo)數(shù)的解析解,再用數(shù)值微分求取原函數(shù)的14 階導(dǎo)數(shù)，并和解析解比較精度。 h=0.05; x=0:h:pi; syms x1; y=sin(x1)/(x12+4*x1+3);%

30、 求各階導(dǎo)數(shù)的解析解與對(duì)照數(shù)據(jù) yy1=diff(y); f1=subs(yy1,x1,x); yy2=diff(yy1); f2=subs(yy2,x1,x); yy3=diff(yy2); f3=subs(yy3,x1,x); yy4=diff(yy3); f4=subs(yy4,x1,x);第51頁(yè)/共103頁(yè) y=sin(x)./(x.2+4*x+3); % 生成已知數(shù)據(jù)點(diǎn) y1,dx1=diff_ctr(y,h,1); subplot(221),plot(x,f1,dx1,y1,:); y2,dx2=diff_ctr(y,h,2); subplot(222),plot(x,f2,dx

31、2,y2,:) y3,dx3=diff_ctr(y,h,3); subplot(223),plot(x,f3,dx3,y3,:); y4,dx4=diff_ctr(y,h,4); subplot(224),plot(x,f4,dx4,y4,:)求最大相對(duì)誤差： norm(y4-f4(4:60)./f4(4:60)ans = 3.5025e-004第52頁(yè)/共103頁(yè)用插值、擬合多項(xiàng)式的求導(dǎo)數(shù) 基本思想：當(dāng)已知函數(shù)在一些離散點(diǎn)上的函數(shù)值時(shí)，該函數(shù)可用插值或擬合多項(xiàng)式來(lái)近似，然后對(duì)多項(xiàng)式進(jìn)行微分求得導(dǎo)數(shù)。 選取x=0附近的少量點(diǎn) 進(jìn)行多項(xiàng)式擬合或插值 g(x)在x=0處的k階導(dǎo)數(shù)為( ,),1,2

32、,1iix yin1121( )nnnng xcxc xc x c()1(0)!0,1,2,knkgckkn 第53頁(yè)/共103頁(yè) 通過(guò)坐標(biāo)變換用上述方法計(jì)算任意x點(diǎn)處的導(dǎo)數(shù)值 令 將g(x)寫(xiě)成z的表達(dá)式 導(dǎo)數(shù)為 可直接用 擬合節(jié)點(diǎn) 得到系數(shù) d=polyfit(x-a,y,length(xd)-1) zxa1121( )( )nnnng xg zdzd zd z d( )( )1( )(0)!0,1,kknkgagdkkn ( )g z(,)iixa yid第54頁(yè)/共103頁(yè) 例：數(shù)據(jù)集合如下： xd: 0 0.2000 0.4000 0.6000 0.8000 1.000 yd: 0.

33、3927 0.5672 0.6982 0.7941 0.8614 0.9053計(jì)算x=a=0.3處的各階導(dǎo)數(shù)。 xd= 0 0.2000 0.4000 0.6000 0.8000 1.000; yd=0.3927 0.5672 0.6982 0.7941 0.8614 0.9053; a=0.3;L=length(xd); d=polyfit(xd-a,yd,L-1);fact=1; for k=1:L-1;fact=factorial(k),fact;end deriv=d.*factderiv = 1.8750 -1.3750 1.0406 -0.9710 0.6533 0.6376第55

34、頁(yè)/共103頁(yè) 建立用擬合（插值）多項(xiàng)式計(jì)算各階導(dǎo)數(shù)的poly_drv.mfunction der=poly_drv(xd,yd,a)m=length(xd)-1;d=polyfit(xd-a,yd,m);c=d(m:-1:1); 去掉常數(shù)項(xiàng)fact(1)=1;for i=2:m; fact(i)=i*fact(i-1);endder=c.*fact; 例： xd= 0 0.2000 0.4000 0.6000 0.8000 1.000; yd=0.3927 0.5672 0.6982 0.7941 0.8614 0.9053; a=0.3; der=poly_drv(xd,yd,a)der

35、= 0.6533 -0.9710 1.0406 -1.3750 1.8750第56頁(yè)/共103頁(yè)二元函數(shù)的梯度計(jì)算 格式： 若z矩陣是建立在等間距的形式生成的網(wǎng)格基礎(chǔ)上，則實(shí)際梯度為,( )xyffgradient z/,/xxyyffxffy( ,)zfx y第57頁(yè)/共103頁(yè) 例：計(jì)算梯度，繪制引力線圖： x,y=meshgrid(-3:.2:3,-2:.2:2); z=(x.2-2*x).*exp(-x.2-y.2-x.*y); fx,fy=gradient(z); fx=fx/0.2; fy=fy/0.2; contour(x,y,z,30); hold on; quiver(x,y

36、,fx,fy)%繪制等高線與引力線圖第58頁(yè)/共103頁(yè) 繪制誤差曲面: zx=-exp(-x.2-y.2-x.*y).*(-2*x+2+2*x.3+x.2.*y-4*x.2-2*x.*y); zy=-x.*(x-2).*(2*y+x).*exp(-x.2-y.2-x.*y); surf(x,y,abs(fx-zx); axis(-3 3 -2 2 0,0.08) figure; surf(x,y,abs(fy-zy); axis(-3 3 -2 2 0,0.11)建立一個(gè)新圖形窗口第59頁(yè)/共103頁(yè) 為減少誤差，對(duì)網(wǎng)格加密一倍： x,y=meshgrid(-3:.1:3,-2:.1:2);

37、 z=(x.2-2*x).*exp(-x.2-y.2-x.*y); fx,fy=gradient(z); fx=fx/0.1; fy=fy/0.1; zx=-exp(-x.2-y.2-x.*y).*(-2*x+2+2*x.3+x.2.*y-4*x.2-2*x.*y); zy=-x.*(x-2).*(2*y+x).*exp(-x.2-y.2-x.*y); surf(x,y,abs(fx-zx); axis(-3 3 -2 2 0,0.02) figure; surf(x,y,abs(fy-zy); axis(-3 3 -2 2 0,0.06)第60頁(yè)/共103頁(yè)3.4 3.4 數(shù)值積分問(wèn)題 由給

38、定數(shù)據(jù)進(jìn)行梯形求積Sum(2*y(1:end-1,:)+diff(y).*diff(x)/2第61頁(yè)/共103頁(yè) 格式： S=trapz(x,y) 例： x1=0:pi/30:pi; y=sin(x1) cos(x1) sin(x1/2); x=x1 x1 x1; S=sum(2*y(1:end-1,:)+diff(y).*diff(x)/2S = 1.9982 0.0000 1.9995 S1=trapz(x1,y) % 得出和上述完全一致的結(jié)果S1 = 1.9982 0.0000 1.9995第62頁(yè)/共103頁(yè) 例：畫(huà)圖 x=0:0.01:3*pi/2, 3*pi/2; % 這樣賦值能確

39、保 3*pi/2點(diǎn)被包含在內(nèi) y=cos(15*x); plot(x,y)% 求取理論值 syms x, A=int(cos(15*x),0,3*pi/2)A =1/15第63頁(yè)/共103頁(yè)隨著步距h的減小，計(jì)算精度逐漸增加： h0=0.1,0.01,0.001,0.0001,0.00001,0.000001; v=; for h=h0, x=0:h:3*pi/2, 3*pi/2; y=cos(15*x); I=trapz(x,y); v=v; h, I, 1/15-I ;end vv = 0.1000 0.0539 0.0128 0.0100 0.0665 0.0001 0.0010 0.0

40、667 0.0000 0.0001 0.0667 0.0000 0.0000 0.0667 0.0000 0.0000 0.0667 0.0000 format long，v第64頁(yè)/共103頁(yè)單變量數(shù)值積分問(wèn)題求解 梯形公式 格式：(變步長(zhǎng)） y=quad(Fun,a,b) y=quadl(Fun,a,b) % 求定積分 y=quad(Fun,a,b, ) y=quadl(Fun,a,b, ) %限定精度的定積分求解，默認(rèn)精度為106。后面函數(shù)算法更精，精度更高。第65頁(yè)/共103頁(yè) 例：第三種：匿名函數(shù)(MATLAB 7.0)第二種：inline 函數(shù)第一種，一般函數(shù)方法第66頁(yè)/共103

41、頁(yè)函數(shù)定義被積函數(shù)： y=quad(c3ffun,0,1.5)y = 0.9661 用 inline 函數(shù)定義被積函數(shù)： f=inline(2/sqrt(pi)*exp(-x.2),x); y=quad(f,0,1.5)y = 0.9661 運(yùn)用符號(hào)工具箱： syms x, y0=vpa(int(2/sqrt(pi)*exp(-x2),0,1.5),60) y0 = y=quad(f,0,1.5,1e-20) % 設(shè)置高精度，但該方法失效第67頁(yè)/共103頁(yè) 例：提高求解精度： y=quadl(f,0,1.5,1e-20)y = 0.9661 abs(y-y0)ans = .640252284

42、8913892e-16 format long 16位精度 y=quadl(f,0,1.5,1e-20)y =第68頁(yè)/共103頁(yè) 例：求解繪制函數(shù): x=0:0.01:2, 2+eps:0.01:4,4; y=exp(x.2).*(x2); y(end)=0; x=eps, x; y=0,y; fill(x,y,g)為減少視覺(jué)上的誤差，對(duì)端點(diǎn)與間斷點(diǎn)（有跳躍）進(jìn)行處理。第69頁(yè)/共103頁(yè) 調(diào)用quad( ): f=inline(exp(x.2).*(x2)./(4-sin(16*pi*x),x); I1=quad(f,0,4)I1 = 57.76435412500863 調(diào)用quadl(

43、): I2=quadl(f,0,4)I2 = 57.76445016946768 syms x; I=vpa(int(exp(x2),0,2)+int(80/(4-sin(16*pi*x),2,4) I =第70頁(yè)/共103頁(yè)求積公式 為使求積公式得到較高的代數(shù)精度 對(duì)求積區(qū)間a,b，通過(guò)變換 有110( )()nkkkfx dxA fx22babaxt110( )()()222222nbkakb ab aa bb ab aa bf x dxftdtA ft1,0.5773503,1.00000000;2,0.7745967,0.555555560.00000000,0.88888889;kk

44、kkkknxAnxAxA 第71頁(yè)/共103頁(yè) 以n=2的高斯公式為例：function g=gauss2(fun,a,b)h=(b-a)/2;c=(a+b)/2;x=h*(-0.7745967)+c, c, h*0.7745967+c;g=h*(0.55555556*(gaussf(x(1)+gaussf(x(3)+0.88888889*gaussf(x(2);function y=gaussf(x)y=cos(x); gauss2(gaussf,0,1)ans = 0.8415第72頁(yè)/共103頁(yè)基于樣條插值的數(shù)值微積分運(yùn)算 基于樣條插值的數(shù)值微分運(yùn)算 格式： Sd=fnder(S,k)該

45、函數(shù)可以求取S的k階導(dǎo)數(shù)。 格式： Sd=fnder(S,k1,kn)可以求取多變量函數(shù)的偏導(dǎo)數(shù)第73頁(yè)/共103頁(yè) 例： syms x; f=(x2-3*x+5)*exp(-5*x)*sin(x); ezplot(diff(f),0,1), hold on x=0:.12:1; y=(x.2-3*x+5).*exp(-5*x).*sin(x); sp1=csapi(x,y);建立三次樣條函數(shù) dsp1=fnder(sp1,1); fnplt(dsp1,-)繪制樣條圖 sp2=spapi(5,x,y);5階次B樣條 dsp2=fnder(sp2,1); fnplt(dsp2,:); axis(

46、0,1,-0.8,5)第74頁(yè)/共103頁(yè) 例：擬合曲面 x0=-3:.3:3; y0=-2:.2:2; x,y=ndgrid(x0,y0); z=(x.2-2*x).*exp(-x.2-y.2-x.*y); sp=spapi(5,5,x0,y0,z); B樣條dspxy=fnder(sp,1,1); fnplt(dspxy)第75頁(yè)/共103頁(yè)理論方法： syms x y; z=(x2-2*x)*exp(-x2-y2-x*y); ezsurf(diff(diff(z,x),y),-3 3,-2 2)對(duì)符號(hào)變量表達(dá)式做三維表面圖第76頁(yè)/共103頁(yè) 基于樣條插值的數(shù)值積分運(yùn)算 格式： f=fn

47、int(S)其中S為樣條函數(shù)。例：考慮 中較稀疏的樣本點(diǎn),用樣條積分的方式求出定積分及積分函數(shù)。 x=0,0.4,1 2,pi; y=sin(x); sp1=csapi(x,y); a=fnint(sp1,1); 建立三次樣條函數(shù) xx=fnval(a,0,pi); xx(2)-xx(1)ans = 2.01910sin xdx第77頁(yè)/共103頁(yè) sp2=spapi(5,x,y); b=fnint(sp2,1); xx=fnval(b,0,pi); xx(2)-xx(1)ans = 1.9999繪制曲線 ezplot(-cos(t)+2,0,pi); hold on不定積分可上下平移 fnp

48、lt(a,-); fnplt(b,:)第78頁(yè)/共103頁(yè)雙重積分問(wèn)題的數(shù)值解 矩形區(qū)域上的二重積分的數(shù)值計(jì)算 格式： 矩形區(qū)域的雙重積分： y=dblquad(Fun,xm,xM,ym,yM) 限定精度的雙重積分： y=dblquad(Fun,xm,xM,ym,yM， )第79頁(yè)/共103頁(yè) 例：求解 f=inline(exp(-x.2/2).*sin(x.2+y),x,y); y=dblquad(f,-2,2,-1,1)y = 1.57449318974494第80頁(yè)/共103頁(yè) 任意區(qū)域上二元函數(shù)的數(shù)值積分 (調(diào)用工具箱N(xiāo)IT），該函數(shù)指定順序先x后y.第81頁(yè)/共103頁(yè)例 fh=in

49、line(sqrt(1-x.2/2),x); % 內(nèi)積分上限 fl=inline(-sqrt(1-x.2/2),x); % 內(nèi)積分下限 f=inline(exp(-x.2/2).*sin(x.2+y),y,x); % 交換順序的被積函數(shù) y=quad2dggen(f,fl,fh,-1/2,1,eps)y =第82頁(yè)/共103頁(yè) 解析解方法： syms x y i1=int(exp(-x2/2)*sin(x2+y), y, -sqrt(1-x2/2), sqrt(1-x2/2); int(i1, x, -1/2, 1)Warning: Explicit integral could not be

50、 found. In D:MATLAB6p5toolboxsymbolicsymint.m at line 58 ans = int(2*exp(-1/2*x2)*sin(x2)*sin(1/2*(4-2*x2)(1/2), x = -1/2 . 1) vpa(ans) ans = 第83頁(yè)/共103頁(yè)222112211sin()yxyIexy dxdy222221sin()xxyIexy dxdy例：計(jì)算單位圓域上的積分： 先把二重積分轉(zhuǎn)化： syms x y i1=int(exp(-x2/2)*sin(x2+y), x, -sqrt(1-y.2), sqrt(1-y.2);Warning:

51、 Explicit integral could not be found. In D:MATLAB6p5toolboxsymbolicsymint.m at line 58第84頁(yè)/共103頁(yè)對(duì)x是不可積的，故調(diào)用解析解方法不會(huì)得出結(jié)果，而數(shù)值解求解不受此影響。 fh=inline(sqrt(1-y.2),y); % 內(nèi)積分上限 fl=inline(-sqrt(1-y.2),y); % 內(nèi)積分下限 f=inline(exp(-x.2/2).*sin(x.2+y),x,y); %交換順序的被積函數(shù) I=quad2dggen(f,fl,fh,-1,1,eps)Integral did not c

52、onverge-singularity likelyI = 0.53686038269795第85頁(yè)/共103頁(yè)三重定積分的數(shù)值求解 格式: I=triplequad(Fun,xm,xM,ym,yM, zm,zM， ,quadl) 其中quadl為具體求解一元積分的數(shù)值函數(shù),也可選用quad或自編積分函數(shù),但調(diào)用格式要與quadl一致。第86頁(yè)/共103頁(yè) 例： x,y,z), 0, 1, 0, pi, 0, pi,1e-7,quadl)ans = 1.7328第87頁(yè)/共103頁(yè)3.5 3.5 曲線積分與曲面積分的計(jì)算曲線積分及MATLAB求解 第一類(lèi)曲線積分 起源于對(duì)不均勻分布的空間曲線總質(zhì)量的求取.設(shè)空間曲線L的密度函數(shù)為f(x,y,z),則其總質(zhì)量 其中s為曲線上某點(diǎn)的弧長(zhǎng),又稱(chēng)這類(lèi)曲線積分為對(duì)弧長(zhǎng)的曲線積分.第88頁(yè)/共103頁(yè) 數(shù)學(xué)表示 若弧長(zhǎng)表示為( ),( ),( )xx tyy tzz t第89頁(yè)/共103頁(yè)例： syms t; syms a positive; x=a