單擊此處編輯母版標題樣式,,單擊此處編輯母版文本樣式,,第二級,,第三級,,第四級,,第五級,,,*,第四講 微積分問題的計算機求解（上）,微積分問題的解析解,,函數(shù)的級數(shù)展開與級數(shù)求和問題求解,,數(shù)值微分,,數(shù)值積分問題,,曲線積分與曲面積分的計算,,,3.1 微積分問題的解析解,3.1.1 極限問題的解析解,單變量函數(shù)的極限,,,,格式1：,L= limit( fun, x, x0),,,,,格式2：,L= limit( fun, x, x0, ‘left’ 或 ‘right’),,例: 試求解極限問題,,>> 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(x^3)-1)/(1-cos(sqrt(x-sin(x)))),x,0,'right'),,ans =,,12,,在(-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:,,divideByZero" to,,suppress this warning.),,>> plot(x,y,'-',[0],,,[12],'o'),,多變量函數(shù)的極限：,,,,,格式：,L,1,=limit(limit(f,x,x,0,),y,y,0,),,或,L,1,=limit(limit(f,y,y,0,), x,x,0,),,,如果x,0,或y,0,不是確定的值，而是另一個變量的函數(shù)，如x->g(y)，則上述的極限求取順序不能交換。,,例：求出二元函數(shù)極限值,,,,,>> syms x y a;,,>> f=exp(-1/(y^2+x^2)) … *sin(x)^2/x^2*(1+1/y^2)^(x+a^2*y^2);,,>> L=limit(limit(f,x,1/sqrt(y)),y,inf),,L =,,exp(a^2),,3.1.2 函數(shù)導數(shù)的解析解,函數(shù)的導數(shù)和高階導數(shù),,,,格式：,y=diff(fun,x),%求導數(shù)(默認為1階),,,y= diff(fun,x,n),%求n階導數(shù),,,例：,,一階導數(shù)：,,>> syms x; f=sin(x)/(x^2+4*x+3);,,>> f1=diff(f); pretty(f1),,cos(x) sin(x) (2 x + 4),,--------------- - -------------------,,2 2 2,,x + 4 x + 3 (x + 4 x + 3),,,原函數(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,,原函數(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,,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),,多元函數(shù)的偏導：,,,,格式：,f=diff(diff(f,x,m),y,n),,或,f=diff(diff(f,y,n),x,m),,,例： 求其偏導數(shù)并用圖表示。,,,>> syms x y z=(x^2-2*x)*exp(-x^2-y^2-x*y);,,>> zx=simple(diff(z,x)),,zx =,,-exp(-x^2-y^2-x*y)*(-2*x+2+2*x^3+x^2*y-4*x^2-2*x*y),,,>> zy=diff(z,y),,zy =,,(x^2-2*x)*(-2*y-x)*exp(-x^2-y^2-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 3 -2 2 -0.7 1.5]),,>> 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) % 繪制引力線,,例,,,,,>> syms x y z; f=sin(x^2*y)*exp(-x^2*y-z^2);,,>> 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)),,多元函數(shù)的Jacobi矩陣:,,,,,,,,,格式：,J=jacobian(Y,X),,其中，X是自變量構(gòu)成的向量，Y是由各個函數(shù)構(gòu)成的向量。,,例：,,試推導其 Jacobi 矩陣,,,>> 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 ],,,,隱函數(shù)的偏導數(shù)：,,,,,,,,,格式：,F=-diff(f,x,j,)/diff(f,x,i,),,,例：,,,,>> syms x y; f=(x^2-2*x)*exp(-x^2-y^2-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),,3.1.3 積分問題的解析解,不定積分的推導：,,,格式：,F=int(fun,x),,,例：,,用diff() 函數(shù)求其一階導數(shù)，再積分，,檢驗是否可以得出一致的結(jié)果。,,>> syms x; y=sin(x)/(x^2+4*x+3); y1=diff(y);,,>> y0=int(y1); pretty(y0) % 對導數(shù)積分,,sin(x) sin(x),,- 1/2 ------ + 1/2 ------,,x + 3 x + 1,,對原函數(shù)求4 階導數(shù)，再對結(jié)果進行4次積分,,>> y4=diff(y,4);,,>> y0=int(int(int(int(y4))));,,>> pretty(simple(y0)),,,sin(x),,------------,,2,,x + 4 x + 3,,,,例:,證明,,,,>> syms a x; f=simple(int(x^3*cos(a*x)^2,x)),,f = 1/16*(4*a^3*x^3*sin(2*a*x)+2*a^4 *x^4+6*a^2*x^2*cos(2*a*x)-6*a*x*sin(2*a*x)-3*cos(2*a*x)-3)/a^4,,>> f1=x^4/8+(x^3/(4*a)-3*x/(8*a^3))*sin(2*a*x)+...,,(3*x^2/(8*a^2)-3/(16*a^4))*cos(2*a*x);,,>> simple(f-f1) % 求兩個結(jié)果的差,,ans =,,-3/16/a^4,,定積分與無窮積分計算:,,,,,格式：,I=int(f,x,a,b),,,,,,格式：,I=int(f,x,a,inf),,,,例：,,,,>> syms x; I1=int(exp(-x^2/2),x,0,1.5) ％無解,,I1 =,,1/2*erf(3/4*2^(1/2))*2^(1/2)*pi^(1/2),,>> vpa(I1,70),,ans =,,,>> I2=int(exp(-x^2/2),x,0,inf),,I2 =,,1/2*2^(1/2)*pi^(1/2),,,,多重積分問題的MATLAB求解,,例：,,,,,,,>>,syms x y z; f0=-4*z*exp(-x^2*y-z^2)*(cos(x^2*y)-10*cos(x^2*y)*y*x^2+...,,4*sin(x^2*y)*x^4*y^2+4*cos(x^2*y)*x^4*y^2-sin(x^2*y));,,>> f1=int(f0,z);f1=int(f1,y);f1=int(f1,x);,,>> f1=simple(int(f1,x)),,f1 =,,exp(-x^2*y-z^2)*sin(x^2*y),,,,>> f2=int(f0,z); f2=int(f2,x); f2=int(f2,x);,,>> f2=simple(int(f2,y)),,f2 =,,2*exp(-x^2*y-z^2)*tan(1/2*x^2*y)/(1+tan(1/2*x^2*y)^2),,>> simple(f1-f2),,ans =,,0,,順序的改變使化簡結(jié)果不同于原函數(shù)，但其誤差為0，表明二者實際完全一致。這是由于積分順序不同，得不出實際的最簡形式。,,例：,,,>> syms x y z,,>> int(int(int(4*x*z*exp(-x^2*y-z^2),x,0,1),y,0,pi),z,0,pi),,ans =,,(Ei(1,4*pi)+log(pi)+eulergamma+2*log(2))*pi^2*hypergeom([1],[2],-pi^2),,Ei(n,z)為指數(shù)積分，無解析解，但可求其數(shù)值解：,,>> vpa(ans,60),,ans =,,,3.2 函數(shù)的級數(shù)展開與 級數(shù)求和問題求解,3.2.1 Taylor 冪級數(shù)展開,,,3.2.2 Fourier 級數(shù)展開,,,3.2.3 級數(shù)求和的計算,,,3.2.1 Taylor 冪級數(shù)展開,3.2.1.1 單變量函數(shù)的 Taylor 冪級數(shù)展開,,,例:,,,,,>> syms x; f=sin(x)/(x^2+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,,>> taylor(f,x,9,2),,ans =,,,,,>> syms a; taylor(f,x,5,a) % 結(jié)果較冗長，顯示從略,,ans =,,sin(a)/(a^2+3+4*a) +(cos(a)-sin(a)/(a^2+3+4*a)*(4+2*a))/(a^2+3+4*a)*(x-a) +(-sin(a)/(a^2+3+4*a)-1/2*sin(a)-(cos(a)*a^2+3*cos(a)+4*cos(a)*a-4*sin(a)-2*sin(a)*a)/(a^2+3+4*a)^2*(4+2*a))/(a^2+3+4*a)*(x-a)^2+…,,例:對y=sinx進行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*x^3+1/120*x^5-1/5040*x^7,,p =,,x-1/6*x^3+1/120*x^5-1/5040*x^7+1/362880*x^9,,p =,,x-1/6*x^3+1/120*x^5-1/5040*x^7+1/362880*x^9-1/39916800*x^11,,p =,,x-1/6*x^3+1/120*x^5-1/5040*x^7+1/362880*x^9-1/39916800*x^11+1/6227020800*x^13,,,p =,,,3.2.1.2 多變量函數(shù)的Taylor 冪級數(shù)展開,多變量函數(shù) 在,,的Taylor冪級數(shù)的展開,,例：,？？？,,,,>> syms x y; f=(x^2-2*x)*exp(-x^2-y^2-x*y);,,>> F=maple('mtaylor',f,'[x,y]',8),,F =,,mtaylor((x^2-2*x)*exp(-x^2-y^2-x*y),[x, y],8),,>> maple(‘readlib(mtaylor)’);％讀庫，把函數(shù)調(diào)入內(nèi)存,,>> F=maple('mtaylor',f,'[x,y]',8),,F =,,-2*x+x^2+2*x^3-x^4-x^5+1/2*x^6+1/3*x^7+2*y*x^2+2*y^2*x-y*x^3-y^2*x^2-2*y*x^4-3*y^2*x^3-2*y^3*x^2-y^4*x+y*x^5+3/2*y^2*x^4+y^3*x^3+1/2*y^4*x^2+y*x^6+2*y^2*x^5+7/3*y^3*x^4+2*y^4*x^3+y^5*x^2+1/3*y^6*x,,>> syms a; F=maple('mtaylor',f,'[x=1,y=a]',3);,,>> F=maple('mtaylor',f,'[x=a]',3),,F =,,(a^2-2*a)*exp(-a^2-y^2-a*y)+((a^2-2*a)*exp(-a^2-y^2-a*y)*(-2*a-y)+(2*a-2)*exp(-a^2-y^2-a*y))*(x-a)+((a^2-2*a)*exp(-a^2-y^2-a*y)*(-1+2*a^2+2*a*y+1/2*y^2)+exp(-a^2-y^2-a*y)+(2*a-2)*exp(-a^2-y^2-a*y)*(-2*a-y))*(x-a)^2,,3.2.2 Fourier 級數(shù)展開,,,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;,,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);,,end,,if a+b, F=subs(F,x,x-L-a); end ％換回變量區(qū)域,,例:,,,>> 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, -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),,例:,,,,>> 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]=fseries(f,x,n), y1=subs(f1,x,xx); plot(xx,y1),,end,,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),,……,,3.2.3 級數(shù)求和的計算,是在符號工具箱中提供的,,例:計算,,,,>> format long; sum(2.^[0:63]) %數(shù)值計算,,ans =,,1.844674407370955e+019,,>> sum(sym(2).^[0:200]) % 或 syms k; symsum(2^k,0,200),,％把2定義為符號量可使計算更精確,,ans =,,,>> syms k; symsum(2^k,0,200),,ans =,,,例,:試求解無窮級數(shù)的和,,,,>> syms n; s=symsum(1/((3*n-2)*(3*n+1)),n,1,inf),,%采用符號運算工具箱,,s =,,1/3,,>> m=1:10000000; s1=sum(1./((3*m-2).*(3*m+1)));%數(shù)值計算方法,雙精度有效位16,“大數(shù)吃小數(shù)”,無法精確,,>> format long; s1 % 以長型方式顯示得出的結(jié)果,,s1 =,,0.33333332222165,,例:求解,,,>> 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)),,%實際應為log((x+1)/x),,例:求,,,>> syms m n; limit(symsum(1/m,m,1,n)-log(n),n,inf),,ans =,,eulergamma,,,,>> vpa(ans, 70) % 顯示 70 位有效數(shù)字,,ans =,,,,,3.3 數(shù)值微分,3.3.1 數(shù)值微分算法,,向前差商公式：,,,向后差商公式,,,,,,,兩種中心公式：,,,,,,,,,,,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*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*Dt^2);L0=3;,,％數(shù)值計算diff(X)表示數(shù)組X相鄰兩數(shù)的差,,case 3,,dy=(-diff(yx1)+7*diff(yx2)-6*diff(yx3)-6*diff(yx4)+...,,7*diff(yx5)-diff(yx6))/(8*Dt^3); L0=5;,,case 4,,dy = (-diff(yx1)+11*diff(yx2)-28*diff(yx3)+28*… diff(yx4)-11*diff(yx5)+diff(yx6))/(6*Dt^4); L0=5;,,end,,dy=dy(L0+1:end-L0); dx=([1:length(dy)]+L0-2-(n>2))*Dt;,,,調(diào)用格式：,,,y為 等距實測數(shù)據(jù)， dy為得出的導數(shù)向量， dx為相應的自變量向量，dy、dx的數(shù)據(jù)比y短 。,,例：,,求導數(shù)的解析解,再用數(shù)值微分求取原函數(shù)的1～4 階導數(shù)，并和解析解比較精度。,,,>> h=0.05; x=0:h:pi;,,>> syms x1; y=sin(x1)/(x1^2+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);,,>> 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]=diff_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,,3.3.3 用插值、擬合多項式的求導數(shù),基本思想：當已知函數(shù)在一些離散點上的函數(shù)值時，該函數(shù)可用插值或擬合多項式來近似，然后對多項式進行微分求得導數(shù)。,,選取x=0附近的少量點,,,進行多項式擬合或插值,,,g(x)在x=0處的k階導數(shù)為,,,通過坐標變換用上述方法計算任意x點處的導數(shù)值,,令,,將g(x)寫成z的表達式,,,導數(shù)為,,,,可直接用 擬合節(jié)點 得到系數(shù),,,d=polyfit(x-a,y,length(xd)-1),,,例：數(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=polyfit(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,,建立用擬合（插值）多項式計算各階導數(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)=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,,3.3.4 二元函數(shù)的梯度計算,,,格式：,,,若z矩陣是建立在等間距的形式生成的網(wǎng)格基礎(chǔ)上，則實際梯度為,,,例：,,計算梯度，繪制引力線圖：,,>> [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),,%繪制等高線與,,引力線圖,,繪制誤差曲面:,,>> 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]),,％建立一個新圖形窗口,,為減少誤差，對網(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=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]),,3.4 數(shù)值積分問題,4.3.1 由給定數(shù)據(jù)進行梯形求積,Sum((2*y(1:end-1,:)+diff(y)).*diff(x))/2,,,格式：,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,,例：,,畫圖,,>> x=[0:0.01:3*pi/2, 3*pi/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,,,隨著步距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,,>> v,,v =,,0.1000 0.0539 0.0128,,0.0100 0.0665 0.0001,,0.0010 0.0667 0.0000,,0.0001 0.0667 0.0000,,0.0000 0.0667 0.0000,,0.0000 0.0667 0.0000,,>> format long，v,,3.4.2 單變量數(shù)值積分問題求解,梯形公式,,,,,,格式：(變步長）,,,y=quad(Fun,a,b),,y=quadl(Fun,a,b),% 求定積分,,,y=quad(Fun,a,b, ),,y=quadl(Fun,a,b, ),%限定精度的定積分求解，默認精度為10,－6,。后面函數(shù)算法更精，精度更高。,,例：,第三種：匿名函數(shù)(MATLAB 7.0),第二種：inline 函數(shù),,第一種，一般函數(shù)方法,,函數(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(-x^2),0,1.5),60),,y0 =,,,>> y=quad(f,0,1.5,1e-20) % 設(shè)置高精度，但該方法失效,,例：,,,提高求解精度：,,>> y=quadl(f,0,1.5,1e-20),,y =,,0.9661,,>> abs(y-y0),,ans =,,.6402522848913892e-16,,>> format long ％16位精度,,>> y=quadl(f,0,1.5,1e-20),,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'),,％為減少視覺上的誤,,差，對端點與間斷點,,（有跳躍）進行處理。,,調(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(x^2),0,2)+int(80/(4-sin(16*pi*x)),2,4)),,I =,,,3.4.3 Gauss求積公式,為使求積公式得到較高的代數(shù)精度,,,,,,,對求積區(qū)間[a,b]，通過變換,,有,,,以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,,3.4.4 基于樣條插值的數(shù)值微積分運算,基于樣條插值的數(shù)值微分運算,,格式：,,,S,d,=fnder(S,k),,該函數(shù)可以求取S的k階導數(shù)。,,,格式：,,,S,d,=fnder(S,[k,1,,…,k,n,]),,可以求取多變量函數(shù)的偏導數(shù),,,例：,,,,>> syms x; f=(x^2-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([0,1,-0.8,5]),,例：,,,,擬合曲面,,>> 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),,理論方法：,,>> syms x y; z=(x^2-2*x)*exp(-x^2-y^2-x*y);,,>> ezsurf(diff(diff(z,x),y),[-3 3],[-2 2]),,％對符號變量表達式做三維表面圖,,基于樣條插值的數(shù)值積分運算,,格式：,,,f=fnint(S),,其中S為樣條函數(shù)。,,,例：考慮,中較稀疏的樣本點,用樣條積分的方式求出定積分及積分函數(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.0191,,>> 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,,％不定積分可上下平移,,>> fnplt(a,'--');,,>> fnplt(b,':'),,3.4.5 雙重積分問題的數(shù)值解,矩形區(qū)域上的二重積分的數(shù)值計算,,,,格式：,,矩形區(qū)域的雙重積分：,,,y=dblquad(Fun,x,m,,x,M,,y,m,,y,M,),,,限定精度的雙重積分：,,,y=dblquad(Fun,x,m,,x,M,,y,m,,y,M，,),,例：,求解,,,>> f=inline('exp(-x.^2/2).*sin(x.^2+y)','x','y');,,>> y=dblquad(f,-2,2,-1,1),,y =,,1.57449318974494,,任意區(qū)域上二元函數(shù)的數(shù)值積分,(調(diào)用工具箱NIT），該函數(shù)指定順序先x后y.,,,,,,,,,例,>> 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 =,,,解析解方法：,,>> syms x y,,>> i1=int(exp(-x^2/2)*sin(x^2+y), y, -sqrt(1-x^2/2), sqrt(1-x^2/2));,,>> int(i1, x, -1/2, 1),,Warning: Explicit integral could not be found.,,> In D:\MATLAB6p5\toolbox\symbolic\@sym\int.m at line 58,,ans =,,int(2*exp(-1/2*x^2)*sin(x^2)*sin(1/2*(4-2*x^2)^(1/2)), x = -1/2 .. 1),,>> vpa(ans),,ans =,,,例：,計算單位圓域上的積分：,,,,,先把二重積分轉(zhuǎn)化：,>> syms x y,,i1=int(exp(-x^2/2)*sin(x^2+y), x, -sqrt(1-y.^2), sqrt(1-y.^2));,,Warning: Explicit integral could not be found.,,> In D:\MATLAB6p5\toolbox\symbolic\@sym\int.m at line 58,,對x是不可積的，故調(diào)用解析解方法不會得出結(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 converge--singularity likely,,I =,,0.53686038269795,,3.4.6 三重定積分的數(shù)值求解,,,,格式:,,,,I=triplequad(Fun,x,m,,x,M,,y,m,,y,M,,z,m,,z,M， ,,@quadl),,,其中@quadl為具體求解一元積分的數(shù)值函數(shù),也可選用@quad或自編積分函數(shù),但調(diào)用格式要與quadl一致。,,例：,,,,>> triplequad(inline('4*x.*z.*exp(-x.*x.*y-z.*z)', …,,'x','y','z'), 0, 1, 0, pi, 0, pi,1e-7,@quadl),,,ans =,,,1.7328,,3.5 曲線積分與曲面積分的計算,3.5.1 曲線積分及MATLAB求解,,第一類曲線積分,,起源于對不均勻分布的空間曲線總質(zhì)量的求取.設(shè)空間曲線,L,的密度函數(shù)為f(x,y,z),則其總質(zhì)量,,,,其中s為曲線上某點的弧長,又稱這類曲線積分為對弧長的曲線積分.,,數(shù)學表示,,若,,,,,弧長表示為,,例：,,,,,>> syms t; syms a positive; x=a*cos(t); y=a*sin(t); z=a*t;,,,>> I=int(z^2/(x^2+y^2)*sqrt(diff(x,t)^2+diff(y,t)^2+ diff(z,t)^2),t,0,2*pi),,I =,,8/3*pi^3*a*2^(1/2),,,,>> pretty(I),,3 1/2,,8/3 pi a 2,,例：,,,>> x=0:.001:1.2; y1=x; y2=x.^2; plot(x,y1,x,y2),,％繪出兩條曲線,,>> syms x; y1=x; y2=x^2; I1=int((x^2+y2^2)*sqrt(1+diff(y2,x)^2),x,0,1);,,>> I2=int((x^2+y1^2)*sqrt(1+diff(y1,x)^2),x,1,0); I=I2+I1,,I =,,-2/3*2^(1/2)+349/768*5^(1/2)+7/512*log(-2+5^(1/2)),,3.5.1.2 第二類曲線積分,又稱對坐標的曲線積分，起源于變力,,沿曲線 移動時作功的研究,,,,,曲線 亦為向量,若曲線可以由參數(shù)方程表示,,,,則兩個向量的點乘可由這兩個向量直接得出.,,例：求曲線積分,,,,>> syms t; syms a positive; x=a*cos(t); y=a*sin(t);,,>> F=[(x+y)/(x^2+y^2),-(x-y)/(x^2+y^2)]; ds=[diff(x,t);diff(y,t)];,,>> I=int(F*ds,t,2*pi,0) % 正向圓周,,I =,,2*pi,,例：,,,,>> syms x; y=x^2; F=[x^2-2*x*y,y^2-2*x*y]; ds=[1; diff(y,x)];,,>> I=int(F*ds,x,-1,1),,,,I =,,,,-14/15,,曲面積分與MATLAB語言求解,3.5.2.1 第一類曲面積分,,,其中 為小區(qū)域的面積，故又稱為對面積的曲面積分。曲面 由 給出，則該積分可轉(zhuǎn)換成x-y平面的二重積分為,,例：,,,,％四個平面，其中三個被積函數(shù)的值為0，只須計算一個即可。,,>> syms x y; syms a positive; z=a-x-y;,,>> I=int(int(x*y*z*sqrt(1+diff(z,x)^2+ diff(z,y)^2),y,0,a-x),x,0,a),,,,I =,,1/120*3^(1/2)*a^5,,若曲面由參數(shù)方程,,,曲面積分,,例：,,,,,>> syms u v; syms a positive;,,>> x=u*cos(v); y=u*sin(v); z=v;f=x^2*y+z*y^2;,,>> E=simple(diff(x,u)^2+diff(y,u)^2+diff(z,u)^2);,,>> F=diff(x,u)*diff(x,v)+diff(y,u)*diff(y,v)+diff(z,u)* diff(z,v);,,>> G=simple(diff(x,v)^2+diff(y,v)^2+diff(z,v)^2);,,>> I=int(int(f*sqrt(E*G-F^2),u,0,a),v,0,2*pi),,I =,,1/4*a*(a^2+1)^(3/2)*pi^2+1/8*log(-a+(a^2+1)^(1/2)) *pi^2-1/8*(a^2+1)^(1/2)*a*pi^2,,3.5.2.2 第二類曲面積分,又稱對坐標的曲面積分,,,,,,可轉(zhuǎn)化成第一類曲面積分,,,若曲面由參數(shù)方程給出,,,例：,,,的上半部，且積分沿橢球面的上面。,,％引入?yún)?shù)方程 x=a*sin(u)*cos(v); y=b*sin(u)*sin(v); z=c*cos(u)， u[0,pi/2], v[0,2*pi].,,>> syms u v; syms a b c positive;,,>> x=a*sin(u)*cos(v); y=b*sin(u)*sin(v); z=c*cos(u);,,>> A=diff(y,u)*diff(z,v)-diff(z,u)*diff(y,v);,,>> I=int(int(x^3*A,u,0,pi/2),v,0,2*pi),,I =,,2/5*pi*a^3*c*b,,