最新中科院-matlab在科學計算中的應(yīng)用3

1、中科院-matlab在科學計算中的應(yīng)用 3 1 第三章 微積分問題的計算機求解 微積分問題的解析解微積分問題的解析解 函數(shù)的級數(shù)展開與級數(shù)求和問題求解函數(shù)的級數(shù)展開與級數(shù)求和問題求解 數(shù)值微分數(shù)值微分 數(shù)值積分問題數(shù)值積分問題 曲線積分與曲面積分的計算曲線積分與曲面積分的計算 中科院-matlab在科學計算中的應(yīng)用 3 2 3.1 微積分問題的解析解 3.1.1 極限問題的解析解 單變量函數(shù)的極限 格式1： L= limit( fun, x, x0) 格式2： L= limit( fun, x, x0, left 或 right) 中科院-matlab在科學計算中的應(yīng)用 3 3 例: 試求解極

2、限問題 syms x a b; f=x*(1+a/x)x*sin(b/x); L=limit(f,x,inf) L = exp(a)*b 例:求解單邊極限問題 syms x; limit(exp(x3)-1)/(1-cos(sqrt(x-sin(x),x,0,right) ans = 12 中科院-matlab在科學計算中的應(yīng)用 3 4 在(-0.1,0.1)區(qū)間繪制出函數(shù)曲線： x=-0.1:0.001:0.1; y=(exp(x.3)-1)./(1-cos(sqrt(x-sin(x); Warning: Divide by zero. (Type warning off MATLAB: d

3、ivideByZero to suppress this warning.)? plot(x,y,-,0, 12,o) 中科院-matlab在科學計算中的應(yīng)用 3 5 多變量函數(shù)的極限： 格式： L1=limit(limit(f,x,x0),y,y0) 或 L1=limit(limit(f,y,y0), x,x0) 如果x0 或y0不是確定的值，而是另一個 變量的函數(shù)，如x-g(y)，則上述的極限求 取順序不能交換。 中科院-matlab在科學計算中的應(yīng)用 3 6 例：求出二元函數(shù)極限值 syms x y a; f=exp(-1/(y2+x2) *sin(x)2/x2*(1+1/y2)(x+a

4、2*y2); L=limit(limit(f,x,1/sqrt(y),y,inf) L = exp(a2) 中科院-matlab在科學計算中的應(yīng)用 3 7 3.1.2 函數(shù)導數(shù)的解析解 函數(shù)的導數(shù)和高階導數(shù) 格式： y=diff(fun,x) %求導數(shù) y= diff(fun,x,n) %求n階導數(shù) 例： 一階導數(shù)： syms x; f=sin(x)/(x2+4*x+3); f1=diff(f); pretty(f1) 中科院-matlab在科學計算中的應(yīng)用 3 8 cos(x) sin(x) (2 x + 4) - - - 2 2 2 x + 4 x + 3 (x + 4 x + 3) 原函

5、數(shù)及一階導數(shù)圖： x1=0:.01:5; y=subs(f, x, x1); y1=subs(f1, x, x1); plot(x1,y,x1,y1,:) 更高階導數(shù)： tic, diff(f,x,100); toc elapsed_time = 4.6860 中科院-matlab在科學計算中的應(yīng)用 3 9 原函數(shù)4階導數(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) (x + 4 x + 3) 3

6、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) 中科院-matlab在科學計算中的應(yīng)用 3 10 多元函數(shù)的偏導： 格式： f=diff(diff(f,x,m),y,n) 或

7、f=diff(diff(f,y,n),x,m) 例： 求其偏導數(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) 中科院-matlab在科學計算中的應(yīng)用 3 11 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); s

8、urf(x,y,z), axis(-3 3 -2 2 -0.7 1.5) 中科院-matlab在科學計算中的應(yīng)用 3 12 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); % 偏導 的數(shù)值解 quiver(x,y,zx,zy) % 繪制引力線 中科院-matlab在科學計算中的應(yīng)用 3 13 例 syms x y z; f=sin(x2*y)*exp(-x2*y-z2

9、); df=diff(diff(diff(f,x,2),y),z); df=simple(df); pretty(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 2 sin(x y) x y+ 4 cos(x y) x y - sin(x y) 中科院-matlab在科學計算中的應(yīng)用 3 14 多元函數(shù)的Jacobi矩陣: 格式：J=jacobian(Y,X) 其中，X是自變量構(gòu)成的向量，Y是由各個函數(shù)構(gòu)成的 向量。 中科院-matlab在科學計算中的應(yīng)用 3 15 例： 試推導其 Ja

10、cobi 矩陣 syms r theta phi; x=r*sin(theta)*cos(phi); y=r*sin(theta)*sin(phi); z=r*cos(theta); 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 中科院-matlab

11、在科學計算中的應(yīng)用 3 16 隱函數(shù)的偏導數(shù)： 格式：F=-diff(f,xj)/diff(f,xi) 中科院-matlab在科學計算中的應(yīng)用 3 17 例： syms x y; f=(x2-2*x)*exp(-x2-y2-x*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) 中科院-matlab在科學計算中的應(yīng)用 3 18 3.1.3 積分問題的解析解 不定積分的推導： 格式： F=int(fun,x) 例： 用diff() 函數(shù)求其一

12、階導數(shù)，再積分，檢驗是否可以 得出一致的結(jié)果。 syms x; y=sin(x)/(x2+4*x+3); y1=diff(y); y0=int(y1); pretty(y0) % 對導數(shù)積分 sin(x) sin(x) - 1/2 - + 1/2 - x + 3 x + 1 中科院-matlab在科學計算中的應(yīng)用 3 19 對原函數(shù)求對原函數(shù)求4 4 階導數(shù)，再對結(jié)果進行階導數(shù)，再對結(jié)果進行4 4次積分次積分 y4=diff(y,4); y0=int(int(int(int(y4); pretty(simple(y0) sin(x) - 2 x + 4 x + 3 中科院-matlab在科學計

13、算中的應(yīng)用 3 20 例:說明明 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)+. (3*x2/(8*a2)-3/(16*a4)*cos(2*a*x); simple(f-f1) % 求兩個結(jié)果的差 ans = -3/16/a4 中科院-matlab在科學計算中的應(yīng)用 3 21 定積分與無窮積分計算

14、: 格式： I=int(f,x,a,b) 格式： I=int(f,x,a,inf) 中科院-matlab在科學計算中的應(yīng)用 3 22 例： syms x; I1=int(exp(-x2/2),x,0,1.5) 無解無解 I1 = 1/2*erf(3/4*2(1/2)*2(1/2)*pi(1/2) vpa(I1,70) ans = 1.6925342362935327229308528 I2=int(exp(-x2/2),x,0,inf) I2 = 1/2*2(1/2)*pi(1/2) 2 / 2 ( ) x fxe 2 0 2 ( ) x t erf xe dt 中科院-matlab在科學計算

15、中的應(yīng)用 3 23 多重積分問題的MATLAB求解 例： syms x y z; f0=-4syms x y z; f0=-4* *z z* *exp(-x2exp(-x2* *y-z2)y-z2)* *(cos(x2(cos(x2* *y)-y)- 1010* *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-sin(x2y2-sin(x2* *y);y); f1=int(f0, f1=int(f0,z z);f

16、1=int(f1,);f1=int(f1,y y);f1=int(f1,);f1=int(f1,x x);); f1=simple(int(f1, f1=simple(int(f1,x x) f1 =f1 = exp(-x2 exp(-x2* *y-z2)y-z2)* *sin(x2sin(x2* *y)y) 中科院-matlab在科學計算中的應(yīng)用 3 24 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) ? f

17、2=sin(x2*y)/exp(y*x2 + z2) simple(f1-f2) ans = 0 順序的改變使化簡結(jié)果不同于原函數(shù)，但 其誤差為0，表明二者實際完全一致。這是由 于積分順序不同，得不出實際的最簡形式。 中科院-matlab在科學計算中的應(yīng)用 3 25 例： 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*hypergeo m(1,2,-pi2) ? Ei(n,z)為指數(shù)積分，無解析解，但可求其數(shù)值

18、解： vpa(ans,60) ans = 3.14127228346146476723483588 中科院-matlab在科學計算中的應(yīng)用 3 26 中科院-matlab在科學計算中的應(yīng)用 3 27 3.2 函數(shù)的級數(shù)展開與 級數(shù)求和問題求解 3.2.1 Taylor 冪級數(shù)展開 3.2.2 Fourier 級數(shù)展開 3.2.3 級數(shù)求和的計算 中科院-matlab在科學計算中的應(yīng)用 3 28 3.2.1 Taylor 冪級數(shù)展開 3.2.1.1 單變量函數(shù)的 Taylor 冪級數(shù)展開 中科院-matlab在科學計算中的應(yīng)用 3 29 中科院-matlab在科學計算中的應(yīng)用 3 30 例: s

19、yms x; f=sin(x)/(x2+4*x+3); y1=taylor(f,x,9); pretty(y1) 23 34 4087 3067 515273 386459 1/3 x - 4/9 x2 + - x3 - - x4 + -x5 - - x6 +- x7 - - x8 54 81 9720 7290 1224720 918540 中科院-matlab在科學計算中的應(yīng)用 3 31 taylor(f,x,9,2) ans = 1/15*sin(2)+(1/15*cos(2)-8/225*sin(2)*(x-2)+ (-127/6750*sin(2)- 8/225*cos(2)*(x-

20、2)2 +(23/6750*cos(2)+628/50625*sin(2)*(x-2)3 +(- 15697/6075000*sin(2)+28/50625*cos(2)*(x-2)4 +(203/6075000*cos(2)+6277/11390625*sin(2)*(x-2)5 +(- 585671/2733750000*sin(2)-623/11390625*cos(2)*(x-2)6 +(262453/19136250000*cos(2)+397361/5125781250*sin(2)*(x-2)7 +(- 875225059/34445250000000*sin(2)-131623

21、/35880468750*cos(2)*(x-2)8 syms a; taylor(f,x,5,a) % 結(jié)果較冗長，顯示從略 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+ 中科院-matlab在科學計算中的應(yīng)用 3 32 例:對y=sin

22、x進行Taylor冪級數(shù)展開,并觀察不同階次的近 似效果。 x0=-2*pi:0.01:2*pi; y0=sin(x0); 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) end p = x-1/6*x3+1/120*x5-1/5040*x7 p = x-1/6*x3+1/120*x5-1/5040*x7+1/362880*x9 p = x-1/6*x3+1/120*x5-1/5040*

23、x7+1/362880*x9-1/39916800*x11 p = x-1/6*x3+1/120*x5-1/5040*x7+1/362880*x9- 1/39916800*x11+1/6227020800*x13 中科院-matlab在科學計算中的應(yīng)用 3 33 p = x-1/6*x3+1/120*x5-1/5040*x7+1/362880*x9- 1/39916800*x11+1/6227020800*x13-1/00*x15 中科院-matlab在科學計算中的應(yīng)用 3 34 3.2.1.2 多變量函數(shù)的Taylor 冪級數(shù)展開 多變量函數(shù) 在 的Taylor冪級數(shù)的展開 12 (,) n

24、 f x xx 12 (,) n a aa 中科院-matlab在科學計算中的應(yīng)用 3 35 例：？ 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) 中科院-matlab在科學計算中的應(yīng)用 3 36 maple(readlib(mtaylor);讀庫，把函數(shù)調(diào)入內(nèi)存 F=maple(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-y

25、2*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*y 2*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)*(-2*a- y)+(2*a-2)*exp(-a2-y2-a*y)*(x-a)+(a2-2*a)*exp(-a2-y2- a*y)*(-1+2

26、*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 中科院-matlab在科學計算中的應(yīng)用 3 37 3.2.2 Fourier 級數(shù)展開 中科院-matlab在科學計算中的應(yīng)用 3 38 function A,B,F=fseries(f,x,n,a,b) if nargin=3, a=-pi; b=pi; end L=(b-a)/2; if a+b, f=subs(f,x,x+L+a); end變量區(qū)域互換 A=int(f,x,-L,L)/L; B=; F=A/2; %計算a0 for i=1:

27、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); end if a+b, F=subs(F,x,x-L-a); end 換回變量區(qū)域 中科院-matlab在科學計算中的應(yīng)用 3 39 例: 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/1

28、6, -12/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) 中科院-matlab在科學計算中的應(yīng)用 3 40 例: syms x; f=abs(x)/x; % 定義方波信號 xx=-pi:pi/200:pi; xx=xx(xx=0); xx=sort(xx,- eps,eps); % 剔除零點 yy=subs(f,x,xx); plot(xx,yy,r-.), hold on % 繪 制出理論值并保持坐標系 for n=2:20 a,b,f1=fser

29、ies(f,x,n), y1=subs(f1,x,xx); plot(xx,y1) end 中科院-matlab在科學計算中的應(yīng)用 3 41 a = 0, 0, 0 b = 4/pi, 0 f1 = 4/pi*sin(x) a = 0, 0, 0, 0 b = 4/pi, 0, 4/3/pi f1 = 4/pi*sin(x)+4/3/pi*sin(3*x) 中科院-matlab在科學計算中的應(yīng)用 3 42 3.2.3 級數(shù)求和的計算 是在符號工具箱中提供的 中科院-matlab在科學計算中的應(yīng)用 3 43 例:計算 format long; sum(2.0:63) %數(shù)值計算 ans = 1.

30、844674407370955e+019 sum(sym(2).0:200) % 或 syms k; symsum(2k,0,200) 把2定義為符號量可使計算更精確更精確 ans = 3298846823252565585670602751 syms k; symsum(2k,0,200) ans = 3298846823252565585670602751 中科院-matlab在科學計算中的應(yīng)用 3 44 例:試求解無窮級數(shù)的和 syms n; s=symsum(1/(3*n-2)*(3*n+1),n,1,inf) %采用符號運算工具箱 s = 1/3 m=1:10000000; s1=s

31、um(1./(3*m-2).*(3*m+1);%數(shù) 值計算方法,雙精度有效位16,“大數(shù)吃小數(shù)”,無法精 確 format long; s1 % 以長型方式顯示得出的結(jié)果 s1 = 0.33333332222165 中科院-matlab在科學計算中的應(yīng)用 3 45 例:求解 syms n x s1=symsum(2/(2*n+1)*(2*x+1)(2*n+1),n, 0,inf); simple(s1) % 對結(jié)果進行化簡，MATLAB 6.5 及以前版本因本身 bug 化簡很麻煩 ans = log(2*x+1)2)(1/2)+1)/(2*x+1)2)(1/2)-1) %實際應(yīng)為log(x+

32、1)/x) 中科院-matlab在科學計算中的應(yīng)用 3 46 例:求 syms m n; limit(symsum(1/m,m,1,n)-log(n),n,inf) ans = eulergamma vpa(ans, 70) % 顯示 70 位有效數(shù)字 ans = .577286824888677 中科院-matlab在科學計算中的應(yīng)用 3 47 符號函數(shù)計算器 單變量符號函數(shù)計算器 Taylor 逼近計算器 中科院-matlab在科學計算中的應(yīng)用 3 48 單變量符號函數(shù)計算器（1/3） 在命令窗口中執(zhí)行 funtool 即可調(diào)出單變 量符號函數(shù)計算器。單變量符號函數(shù)計 算器用于對單變量函數(shù)

33、進行操作，可以 對符號函數(shù)進行化簡、求導、繪制圖形 等。該工具的界面如圖所示。 函數(shù) f 的 圖形窗口 函數(shù) g 的 圖形窗口 控制窗口 中科院-matlab在科學計算中的應(yīng)用 3 49 單變量符號函數(shù)計算器（2/3） 1輸入框的功能 如圖： 函數(shù)函數(shù) f 的編輯框的編輯框 函數(shù)函數(shù) g 的編輯框的編輯框 顯示繪制顯示繪制 f 和和 g 的圖像的的圖像的 x 區(qū)間區(qū)間 用于修改用于修改 f 的常數(shù)因子的常數(shù)因子 0 函數(shù) f 自 身的操作 函數(shù) f 與常 數(shù) a 的操作 函數(shù) f 與函 數(shù) g 的操作 系統(tǒng)操作 中科院-matlab在科學計算中的應(yīng)用 3 50 單變量符號函數(shù)計算器（3/3）

34、單變量符號函數(shù)計算器應(yīng)用示例 在 f 函數(shù)輸入欄 中輸入 cos(x3) cos(x3) (1+x2) 在 g 函數(shù)輸入欄 中輸入 (1+x2) 點擊 中科院-matlab在科學計算中的應(yīng)用 3 51 Taylor 逼近計算器 Taylor 逼近計算器用于實現(xiàn)函數(shù)的 taylor 逼近。在命令窗口中輸入 taylortool，調(diào) 出Taylor 逼近計算器，界面及功能如圖。 輸入待逼 近的函數(shù) 輸入擬合函 數(shù)的階數(shù) 級數(shù)的展開點， 默認為 0 輸入擬 合區(qū)間 中科院-matlab在科學計算中的應(yīng)用 3 52 MAPLE 函數(shù)的調(diào)用 maple 函數(shù)的使用 mfun 函數(shù)的使用 中科院-matl

35、ab在科學計算中的應(yīng)用 3 53 maple 函數(shù)的使用 maple 是符號工具箱中的一個通用命令，使用它可以實 現(xiàn)對 MAPLE 中大部分函數(shù)的調(diào)用。其使用格式為： r = maple(statement)，其中 statement 為符合 MAPLE 語法的可 執(zhí)行語句的字符串，該命令將 statement 傳遞給 MAPLE，該 命令的輸出結(jié)果也符合 MAPLE 的語法； r = maple(function,arg1,arg2,.)，該函數(shù)調(diào)用引號中的函數(shù)， 并接受指定的參數(shù)，相當于 MAPLE 語句 function(arg1,arg2,.)； r, status = maple(.

36、)，返回函數(shù)的運行狀態(tài)，如果函數(shù)運行成 功，則 status 為 0，r 為運行結(jié)果；如果函數(shù)運行失敗，則 status 為一個正數(shù)，r 為相應(yīng)的錯誤信息； maple(traceon) 或者 maple trace on，輸出 MAPLE 函數(shù)運行中 的所有子表達式和運行結(jié)果； maple(traceoff) 或 maple trace off，不顯示中間過程。 中科院-matlab在科學計算中的應(yīng)用 3 54 mfun 函數(shù)的使用 mfun 函數(shù)用于對 maple 函數(shù)進行數(shù)字評估。該函數(shù)的 調(diào)用格式為： Y = mfun(function,par1,par2,par3,par4)。 該語

37、句對指定的數(shù)學函數(shù)進行評估。其中 function 指 定待評估的函數(shù)，par1、par2 等為 function 的參數(shù)， 為待評估的數(shù)值，其維數(shù)有 function 函數(shù)的參數(shù)類型 確定。在該語句中最多可以設(shè)置四個參數(shù)，最后一個 參數(shù)可以為矩陣。 用戶可以通過 help mfunlist 查看 MATLAB 中 mfun 可 以調(diào)用的函數(shù)列表，另外，可以通過 mhelp function 查 看指定函數(shù)的相關(guān)信息。 中科院-matlab在科學計算中的應(yīng)用 3 55 3.3 數(shù)值微分 0 ()( ) ( )lim h f xhf x fx h 差商型求導公式 由導數(shù)定義 1 ()( ) (

38、) 2 ( )() ( ) 3 ()() ( ) 2 fxhfx fx h fxfxh fx h fxhfxh fx h （ ）向前差商公式 （ ）向后差商公式 （ ）中心差商公式 （中點方法 ） x-h x x+h B C A T f(x) 中科院-matlab在科學計算中的應(yīng)用 3 56 3.3.1 數(shù)值微分算法 向前差商公式： 向后差商公式 中科院-matlab在科學計算中的應(yīng)用 3 57 兩種中心公式： 234 234 2 ( )( )( )/ 2!( )/3!() ( ) 2 ( )( )( )/ 2!( )/3!() 2 ( )( ) 3! f xtfxt fxt fot f x

39、t f xtfxt fxt fot t t fxf % 中科院-matlab在科學計算中的應(yīng)用 3 58 中科院-matlab在科學計算中的應(yīng)用 3 59 中科院-matlab在科學計算中的應(yīng)用 3 60 3.3.2 中心差分方法及其 MATLAB 實現(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 = (-diff(yx1)+7*d

40、iff(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ù)值計算diff(X)表示數(shù)組X相鄰兩數(shù)的差 中科院-matlab在科學計算中的應(yīng)用 3 61 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(y

41、x2)-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ù)據(jù)， dy為得出的導數(shù)向量， dx為 相應(yīng)的自變量向量，dy、dx的數(shù)據(jù)比y短 。 ,_( , ) yx dddiffctr yt n 中科院-matlab在科學計算中的應(yīng)用 3 62 例： 求導數(shù)的解析解,再用數(shù)值微分求取原函數(shù)的1 4 階導數(shù)，并和解析解比較精度。 h=0.05; x=0:h:pi; syms x1;

42、y=sin(x1)/(x12+4*x1+3); % 求各階導數(shù)的解析解與對照數(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); 中科院-matlab在科學計算中的應(yīng)用 3 63 y=sin(x)./(x.2+4*x+3); % 生成已知數(shù)據(jù)點 y1,dx1=diff_ctr(y,h,1); subplot(221),plot(x,f1,dx1,y1,:); y2,dx2=di

43、ff_ctr(y,h,2); subplot(222),plot(x,f2,dx2,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,:) 求最大相對誤差： norm(y4- f4(4:60)./f4(4:60) ans = 3.5025e-004 中科院-matlab在科學計算中的應(yīng)用 3 64 3.3.3 用插值、擬合多項式的求導數(shù) 基本思想：當已知函數(shù)在一些離散點上的函 數(shù)值時，該函數(shù)可用插值或擬合多項式

44、來近 似，然后對多項式進行微分求得導數(shù)。 選取x=0附近的少量點 進行多項式擬合或插值 g(x)在x=0處的k階導數(shù)為 ( ,),1,2,1 ii x yin 1 121 ( ) nn nn g xcxc xc x c () 1 (0)!0,1,2, k nk gckkn 中科院-matlab在科學計算中的應(yīng)用 3 65 通過坐標變換用上述方法計算任意x點處的導 數(shù)值 令 將g(x)寫成z的表達式 導數(shù)為 可直接用 擬合節(jié)點 得到系數(shù) d=polyfit(x-a,y,length(xd)-1) zxa 1 121 ( )( ) nn nn g xg zdzd zd z d ( )( ) 1 (

45、 )(0)!0,1, kk nk gagdkkn ( )g z(,) ii xa y i d 中科院-matlab在科學計算中的應(yīng)用 3 66 例：數(shù)據(jù)集合如下： 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 計算x=a=0.3處的各階導數(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=polyfi

46、t(xd-a,yd,L-1);fact=1; for k=1:L-1;fact=factorial(k),fact;end deriv=d.*fact deriv = 1.8750 -1.3750 1.0406 -0.9710 0.6533 0.6376 中科院-matlab在科學計算中的應(yīng)用 3 67 建立用擬合（插值）多項式計算各階導數(shù)的 poly_drv.m function der=poly_drv(xd,yd,a) m=length(xd)-1; d=polyfit(xd-a,yd,m); c=d(m:-1:1); 去掉常數(shù)項 fact(1)=1;for i=2:m; fact(i)

47、=i*fact(i-1);end der=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 = 0.6533 -0.9710 1.0406 -1.3750 1.8750 中科院-matlab在科學計算中的應(yīng)用 3 68 3.3.4 二元函數(shù)的梯度計算 格式： 若z矩陣是建立在等間距的形式生成的網(wǎng)格基 礎(chǔ)上，則實際梯度為 ,( ) xy ffgradient z /,/ xxyy ff

48、xffy ( ,)zfx y 中科院-matlab在科學計算中的應(yīng)用 3 69 例： 計算梯度，繪制引力線圖： 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,fx,fy) %繪制等高線與 引力線圖 中科院-matlab在科學計算中的應(yīng)用 3 70 繪制誤差曲面: zx=-exp(-x.2-y.2-x.*y).*(-2*x+2+2*x.3+x.2.*y-4

49、*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) 建立一個新圖形窗口 中科院-matlab在科學計算中的應(yīng)用 3 71 為減少誤差，對網(wǎng)格加密一倍： x,y=meshgrid(-3:.1:3,-2:.1:2); z=(x.2-2*x).*exp(-x.2-y.2-x.*y); fx,fy=gradient(z); fx=fx/0.1; fy

50、=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) 中科院-matlab在科學計算中的應(yīng)用 3 72 3.4 數(shù)值積分問題 4.3.1 由給定數(shù)據(jù)進行梯形求積 中科院-matlab在科學計算中的應(yīng)用 3 73 Sum(2*y(1:end

51、-1,:)+diff(y).*diff(x)/2 中科院-matlab在科學計算中的應(yīng)用 3 74 格式： 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)/2 S = 1.9982 0.0000 1.9995 S1=trapz(x1,y) % 得出和上述完全一致的結(jié)果 S1 = 1.9982 0.0000 1.9995 中科院-matlab在科學計算中的應(yīng)用 3 75 例： 畫圖 x=0:0.01:3*pi/2, 3*pi/

52、2; % 這樣賦值能確保 3*pi/2點 被包含在內(nèi) y=cos(15*x); plot(x,y) % 求取理論值 syms x, A=int(cos (15*x),0,3*pi/2) A = 1/15 中科院-matlab在科學計算中的應(yīng)用 3 76 隨著步距h的減小，計算精度逐漸增加： 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 format long; v v = 0.100

53、 0.076 0.591 0.000 0.584 0.083 0.000 0.004 0.663 0.000 0.667 0.000 0.000 0.167 0.500 0.000 0.542 0.125 中科院-matlab在科學計算中的應(yīng)用 3 77 3.4.2 單變量數(shù)值積分問題求解 梯形公式 格式：(變步長）(Fun:函數(shù)的字符串變量） y=quad(Fun,a,b) y=quadl(Fun,a,b) % 求定積分 y=quad(Fun,a,b, ) y=quadl(Fun,a,b, ) %限定精 度的定積分求解，默認精度為106。后面函數(shù) 算法更精確，精度更高。 中科院-matlab

54、在科學計算中的應(yīng)用 3 78 例： 第三種：匿名函數(shù)(MATLAB 7.0) 第二種：inline 函數(shù) 第一種，一般函數(shù)方法 中科院-matlab在科學計算中的應(yīng)用 3 79 函數(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 運用符號工具箱： syms x, y0=vpa(int(2/sqrt(pi)*exp(-x2),0,1.5),60) y0 = .966179499943257461473285

55、749 y=quad(f,0,1.5,1e-20) % 設(shè)置高精度，但該方法失效 中科院-matlab在科學計算中的應(yīng)用 3 80 提高求解精度： y=quadl(f,0,1.5,1e-20) y = 0.9661 abs(y-y0) ans = .64892e-16 format long 16位精度 y=quadl(f,0,1.5,1e-20) y = 0.96610514647531 中科院-matlab在科學計算中的應(yīng)用 3 81 例：求解 繪制函數(shù): x=0:0.01:2, 2+eps:0.01:4,4; y=exp(x.2).*(x2); y(end)=0; x=eps, x; y

56、=0,y; fill(x,y,g) 為減少視覺上的誤 差，對端點與間斷點 （有跳躍）進行處理。 中科院-matlab在科學計算中的應(yīng)用 3 82 調(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( ): 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 = 57.76445333 中科院-mat

57、lab在科學計算中的應(yīng)用 3 83 3.4.3 Gauss求積公式 為使求積公式得到較高的代數(shù)精度 對求積區(qū)間a,b，通過變換 有 1 1 0 ( )() n kk k fx dxA fx 22 baba xt 1 1 0 ( )()() 222222 n b k a k b ab aa bb ab aa b f x dxftdtA ft 0101 0202 11 1,0.5773503,1.00000000; 2,0.7745967,0.55555556 0.00000000,0.88888889; nxxAA nxxAA xA 中科院-matlab在科學計算中的應(yīng)用 3 84 以n=2的高

58、斯公式為例： 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 22 baba xt 0 ( )() 222 n b k a k b ab aa b f x dxA ft 021 021 2,0.7745967,0.0

59、0000000 0.55555556,0.88888889; nxxx AAA 中科院-matlab在科學計算中的應(yīng)用 3 85 3.4.4 雙重積分問題的數(shù)值解 矩形區(qū)域上的二重積分的數(shù)值計算 格式： 矩形區(qū)域的雙重積分： y=dblquad(Fun,xm,xM,ym,yM) 限定精度的雙重積分： y=dblquad(Fun,xm,xM,ym,yM， ) 中科院-matlab在科學計算中的應(yīng)用 3 86 例：求解 f=inline(exp(-x.2/2).*sin(x.2+y),x,y); y=dblquad(f,-2,2,-1,1) y = 1.57449318974494 中科院-mat

60、lab在科學計算中的應(yīng)用 3 87 任意區(qū)域上二元函數(shù)的數(shù)值積分 (調(diào)用工具箱 NIT），該函數(shù)指定順序先x后y. 中科院-matlab在科學計算中的應(yīng)用 3 88 例 fh=inline(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 = 0.41192954617630 中科院-matlab在科學計算中的應(yīng)用 3 89 解析解方法： syms