拉普拉斯方程與泊松方程

1、第五章 拉普拉斯方程與泊松方程,5.1 二維拉普拉斯方程5.2 三維拉普拉斯方程5.3 泊松方程與格林函數(shù),先求其本征解。設(shè)函數(shù) u(x,t) 具有變量分離形式： ，則上述方程可以寫為：,5.1 二維拉普拉斯方程,5.1.1 矩形區(qū)域的拉普拉斯方程,其分布如下頁左圖所示，右圖則為偏微分方程工具箱算出的結(jié)果。,用偏微分方程工具箱求解Laplace方程的步驟,

2、1、在matlab命令窗中鍵入：pdetool2、在彈出界面中用第二行的左邊5個(gè)按鈕之一選定求解邊界.3、用第二行的左邊第6個(gè)按鈕設(shè)定邊界條件.4、用第二行的左邊第7個(gè)按鈕設(shè)定求解方程 (如選橢圓形并設(shè)定系數(shù)）.5、用第二行的左邊第8、9個(gè)按鈕剖分求解區(qū)域網(wǎng)格.6、用第二行的左邊第11個(gè)按鈕畫圖.,5.1.2 圓形區(qū)域拉普拉斯方程，陽光照射的圓柱,其中：,半徑為a,表面熏黑的長(zhǎng)圓柱體，在溫度為零度的空氣中受到

3、垂直于柱軸的陽光照射，熱流的強(qiáng)度為q, 求柱內(nèi)溫度分布。 因溫度分布式穩(wěn)定的，該問題的定界問題為：,該問題的解析解為 ：,%ex401(p91) clear;a=1; H=1.5; k=0.3; h=0.2; q=5; N=20;r=0:0.05:1; phi=0:pi/30:2*pi; [TH,R]=meshgrid(phi,r); %構(gòu)造網(wǎng)格[X,Y]=pol2cart(TH,R); u=q/H/pi+1/(k

4、+H*a)*q/2*R.*sin(TH);for n=1:N dd=2*q/pi/a^(2*n-1)/(2*n*k+a*H)/(1-4*n^2);zz=dd*R.^(2*n).*cos(2*n*TH); u=u+zz;end;figure(1);surfc(X,Y,u);figure(2);contour3(X,Y,u,20);,求上述溫度分布的程序如下，相應(yīng)的分度分布如下頁上面兩圖，用偏微分方程工具箱所得分布

5、可作為對(duì)比,5.1.3 云與大地之間的電纜,帶電的云與大地之間存在一個(gè)均勻的電場(chǎng)，平行與大地的電纜相當(dāng)于一根無窮長(zhǎng)導(dǎo)體。在平行于電場(chǎng)的方向作垂直于電纜的截面，研究該截面上的電勢(shì)分布。 該問題可以用如下方程加以描述：,選擇偏微分方程工具箱求解，注意：求解區(qū)域是兩區(qū)域之差；區(qū)域內(nèi)選擇橢圓型方程。,5.2 三維拉普拉斯方程,5.2.1 靜電場(chǎng)中的介質(zhì)球的電場(chǎng)分布,在場(chǎng)強(qiáng)為E的均勻靜電場(chǎng)中放置半徑為R0的均勻介質(zhì)球，球的介電常

6、數(shù)(電容率)為ε，求介質(zhì)球內(nèi)外的電場(chǎng)強(qiáng)度分布。 定界問題可以用下列方程加以描述：,圓外ε1=2 ，圓內(nèi)ε2=2,5.2.2 帶有電荷的細(xì)圓環(huán)的電勢(shì)分布,半徑為a 的均勻帶電細(xì)圓環(huán)，電荷線密度為4πε0 /a,取無窮遠(yuǎn)處電勢(shì)為0,求空間任一點(diǎn)的電勢(shì).直接積分法,由于軸對(duì)稱性，可取x=0，可用trapz求先出電勢(shì)，再用gradient(梯度)求出電場(chǎng)，最后用streamline畫出電力線,%ex403(p97) clear;

7、 a=1; b=0.11; y=-4:b:4; z=y; phi=pi*(0:1/100:2); [Y,Z,PHI]=meshgrid(y,z,phi); r=sqrt((0-a*cos(PHI)).^2+(Y-a*sin(PHI)).^2+Z.^2);dV=1./r; V=trapz(dV,3); [Ey,Ez]=-gradient(V,0.5);figure(2);subplot(2,2,1); contour(Y(:,

8、:,1),Z(:,:,1),V,10); subplot(2,2,3); [Sy,Sz]=meshgrid(-4:.2:4,[-0.1,.1]); box on; streamline(Y(:,:,1),Z(:,:,1),Ey,Ez,Sy,Sz);x=0:b:3; th=pi*(0:1/20:2); [X,Y,Z,TH]=ndgrid(x,x,x,th);r=sqrt((X-a*cos(TH)).

9、^2+(Y-a*sin(TH)).^2+Z.^2);dV=1./r; V=trapz(dV,4); [Ex,Ey,Ez]=gradient(-V,0.5); [X,Y,Z]=meshgrid(x); [Sx,Sy,Sz]=meshgrid(0:.5:3,0:.5:3,0.1); subplot(2,2,2); streamslice(X,Y,Z,Ex,Ey,Ez,[],[],0.1); box on; axis([0,3,0,

10、3]);xlabel('Z_0=0.1');x=cos(th);y=sin(th);z=zeros(1,length(th));subplot(2,2,4); plot3(x,y,z,'linewidth',3,'color','r');hold on; axis([-3,3,-3,3,-3,3]);h1=streamline(X,Y,Z,Ex,Ey,Ez,Sx,Sy,

11、Sz);h2=copyobj(h1,gca); rotate(h2,[1,0,0],180,[0 0 0]); h3=copyobj(allchild(gca),gca); rotate(h3,[0,1,0],180,[0 0 0]);,2. 解析解的可視化,該定解問題的解析解----電勢(shì)分布為：,由于對(duì)稱性， 解與角φ無關(guān)，該解析解可以分為圓內(nèi)、圓外兩部分，兩部分分別計(jì)算，再將兩部分相加。程序如下頁；相應(yīng)的等位線在后頁。,%ex

12、402(p94) clear; a=0.5; q=1; x=[-3*a:0.02:3*a]; [X,Y]=meshgrid(x); [theta,r]=cart2pol(Y,X);rout=r; rout(find(routa))=NaN;Uin=q/a; Uout=q./rout; rin=rin/a; rout=a./rout;for k=1:20 fun=legendre(2*k,cos(th

13、eta)); rfun=q/a*squeeze(fun(1,:,:)); ck=(-1)^k*prod(1:2*k)/2^(2*k)/(prod(1:k))^2; ukin=ck*rin.^(2*k); Uin=Uin+ukin.*rfun; ukout=ck*rout.^(2*k+1); Uout=Uout+ukout.*rfun;end;figure(1);contour(X,Y

14、,Uout,20,'r');hold on;contour(X,Y,Uin,15,'b');title('等勢(shì)線');figure(2);surf(X,Y,Uout);hold on;surf(X,Y,Uin);,3. 用偏微分方程工具箱求解,5.2.3 均勻圓盤的引力勢(shì),均勻圓盤的半徑為a，質(zhì)量為M，求它在周圍空間中的引力勢(shì)。定解問題：其中：下面分別用兩種方法求解

15、：解析法繪圖， 2. PDE工具箱求解,泊松方程,1. 解析解的可視化問題的解析解為：考慮到對(duì)稱性，解析解與φ無關(guān)，因此在二維平面上畫出等勢(shì)線，由于解為分段函數(shù)，因此分別算出并用不同顏色作圖。,%ex404; (p100) 圓盤引力勢(shì)；clear; a=0.35; GM=1/4; ri=0:1/200:a; ro=a:1/200:1;th=(0:0.01:2)*pi; z=cos(th); ui=0;

16、 uo=0;for k=0:2:20 fun=legendre(k,z); f=fun(1,:); f0=f(1,51); Ri=2/a^2*(1/(k-1)+1/(k+2)).*ri-2/a/(k-1)*(ri/a).^k; Ro=1/a/(k/2+1)*(a./ro).^(k+1); [R,PH]=meshgrid(Ri,f); ui =ui-GM*f0*R.*PH; [R,PH]=

17、meshgrid(Ro,f); uo=uo-GM*f0*R.*PH; end;[ro,T]=meshgrid(ro,th);Yo=ro.*cos(T);Xo=ro.*sin(T);[ri,T]=meshgrid(ri,th); Yi=ri.*cos(T); Xi=ri.*sin(T);figure(1);contour(Xo,Yo,uo,'b');hold on; axis equal;contour(Xi

18、,Yi,ui,'r'); title('圓盤引力勢(shì)等勢(shì)線');,2. 用PDE工具箱求解 選一長(zhǎng)為3，寬為2的矩形，內(nèi)畫一小矩形表示板，長(zhǎng)(直徑)取為0.8, 厚度取為0.1，中心位于原點(diǎn)。 外邊界條件是Dirichlet型：h=1,r=0; 內(nèi)邊界條件自動(dòng)銜接。 方程選為橢圓型，在圓板內(nèi)參數(shù)取為：c=1,a=0, f=1/0.35^2./sqrt(x.^2+y.

19、^2); 圓板外參數(shù)取為：c=1,a=0, f=0。,5.2.4 環(huán)形電流的磁感應(yīng)強(qiáng)度,半徑為a，通有電流I的圓環(huán)，求其在空間任一點(diǎn)所產(chǎn)生的磁感應(yīng)強(qiáng)度B.為了畫出磁力線，需用到畫場(chǎng)線的專用命令(streamline)。解析解,%ex406; (p104) % 環(huán)形電流的磁感應(yīng)強(qiáng)度clear;a=0.35; R=1;[X,Y]=meshgrid(0:0.1:R);[t,r]=cart2pol(Y,X);ri=r; r

20、i(find(ri>a))=NaN; ro=r; ro(find(ro<a))=NaN;Bri=0; Bqi=0; Bro=0; Bqo=0; for k=1:2:43 fun10=legendre(k,0); P0=fun10(2); fun=legendre(k,cos(t));P=squeeze(fun(1,:,:)); P1=squeeze(fun(2,:,:));

21、 Cr=ri.^(k-1)/a^k; Bri=Bri-1/2*P0*Cr.*P; Dr=(a./ro).^(k+1)./ro; Bro=Bro-1/2*P0*Dr.*P; Cq=1/k*Cr; Bqi=Bqi-1/2*P0*Cq.*P1Dq=1/(k+1)*Dr; Bqo=Bqo+1/2*P0*Dq.*P1end;,%--------------------

22、-------------------將非數(shù)元素令其為零 L1=isnan(Bri); Bri(find(L1==1))=0; Bqi(find(L1==1))=0; L2=isnan(Bro);Bro(find(L2==1))=0; Bqo(find(L2==1))=0; Br=Bri+Bro; Bq=Bqi+Bqo; %內(nèi)外合并； Bx=Br.*sin(t)+Bq.*cos(t); %換為直角坐

23、標(biāo)分量； By=Br.*cos(t)-Bq.*sin(t); vy=0;vx=[0:0.025:a];[Vx,Vy]=meshgrid(vx,vy); figure(1); hold on;%---------------------------------------畫圖streamline(X,Y,Bx,By,Vx,Vy); streamline(-X,Y,-Bx,By,-Vx,Vy);stre

24、amline(-X,-Y,-Bx,-By,-Vx,-Vy);treamline(X,-Y,Bx,-By,Vx,-Vy); box on; title('電流環(huán)的磁力線');figure(2);subplot(2,2,1);surf(X,Y,Br); title('r分量'); view(-130,20);subplot(2,2,2);contour(X,Y,Br,10); title(&

25、#39;r分量等值線');subplot(2,2,3);surf(X,Y,Bq); title('\theta分量');subplot(2,2,4);contour(X,Y,Bq,10); title('\theta分量等值線');,2. 利用橢圓函數(shù)表示的直接積分 根據(jù)比奧-薩伐爾定律：可積分得：其中E,K為橢圓函數(shù)：,%ex407; (p107) % 環(huán)形電流的

26、磁感應(yīng)強(qiáng)度(橢圓函數(shù)法)；clear; a=0.35; R=1; [X,Y]=meshgrid(0:0.1:R);rr2=X.^2+Y.^2; k2=4*a*X./(a^2+r2+2*a*X);[K,E]=ellipke(k2); D=a^2+r2-2*a*X; D=D+eps*(D==0);aE1=(a^2-r2)./D.*E; aE2=(a^2+r2)./D.*E;Bx=2*Y./(X+eps*(X=

27、=0))./sqrt(a^2+r2+2*a*X).*(aE2-K);By=2./sqrt(a^2+r2+2*a*X).*(aE1+K);vy=0; vx=[0.0:0.025:a]; [Vx,Vy]=meshgrid(vx,vy); figure(1); streamline(X,Y,Bx,By,Vx,Vy); hold on; box on; streamline(-X,Y,-Bx,By,-Vx,Vy); streamlin

28、e(X,-Y,Bx,-By,Vx,-Vy);streamline(-X,-Y,-Bx,-By,-Vx,-Vy); title('電流環(huán)的磁力線');figure(2);subplot(2,2,1);surf(X,Y,Bx);title('Bx分量');view(-130,20);subplot(2,2,2);contour(X,Y,Bx,10);title('Bx分量等值線');

29、subplot(2,2,3);surf(X,Y,By);title('By分量');subplot(2,2,4);contour(X,Y,By,10);title('By分量等值線');,5.2.5 柱體內(nèi)溫度場(chǎng)分布之一(Jν的應(yīng)用),半徑為b高度為h的圓柱體，其側(cè)面絕熱，下底溫度保持0度，上底溫度分布為 。求圓柱體內(nèi)穩(wěn)定的溫度分布。定解問題為：其解析解為：；其中

30、 為2階貝塞爾函數(shù)在1階貝塞爾的第n個(gè)0點(diǎn)xn1處的值；而 為0階貝塞爾函數(shù)在1階貝塞爾的第n個(gè)0點(diǎn)xn1處的值取b =1.5，u0=1，h=0.5, 用下列程序畫圖，結(jié)果如后頁所示。,%ex410; (p115) % 溫度場(chǎng)分布；clear; D=[]; N=50; j=1; yy=inline('besselj(1,xx)','xx');for k=1:N;

31、while(besselj(1,j)*besselj(1,j+1)>0), j=j+1;end; %粗估零點(diǎn);q=fzero(yy,j); D=[D,q]; j=j+1;%%精估一階貝塞爾函數(shù)零點(diǎn)；endb=1.5; u0=1; h=0.5; [rho,z]=meshgrid(0:b/30:b,0:h/20:h); A=u0*b^2/(2*h)*z;for n=1:N;aa=4*u0*besselj(2,D(n)

32、)/(D(n)/b)^2;bb=sinh(D(n)*z./b)/sinh(D(n)*h/b);cc=besselj(0,D(n).*rho/b)/(besselj(0,D(n)))^2;A=A-aa*bb.*cc;end;figure(1);contour(0:b/30:b,0:h/20:h,A,25);xlabel(‘radical \rho'); ylabel('high z'); colorb

33、ar;,5.2.6 柱體內(nèi)溫度場(chǎng)分布之二 (Jν的應(yīng)用),半徑為ρ 0高度為L(zhǎng)的圓柱體，其側(cè)面和下底溫度保持0度，上底溫度分布為 。求圓柱體內(nèi)穩(wěn)定的溫度分布。定解問題為：其解析解為：其中 為零階貝塞爾函數(shù)的第n個(gè)零點(diǎn)，J0和J1為0、2階貝塞爾函數(shù)取ρ0 =0.2，L=0.5, 用下列程序畫圖，結(jié)果如后頁所示。,%ex411; (p117) % 溫度場(chǎng)分布之二；clear; D=[];

34、N=50; j=1; yy=inline('besselj(0,xx)','xx');for k=1:N; while(besselj(0,j)*besselj(0,j+1)>0), j=j+1;end;%零點(diǎn)粗估；q=fzero(yy,j); D=[D,q]; j=j+1; %精估零階貝塞爾函數(shù)零點(diǎn)；endr0=0.2; L=0.8; A=0;[r,z]=meshgri

35、d(-r0:0.01:r0,0:0.01:L);for n=1:N;aa=2*r0^2/D(n)/besselj(1,D(n))*(1-4/(D(n))^2);bb=sinh(D(n)*z/r0)./sinh(D(n)*L/r0);cc=besselj(0,D(n)*r/r0);A=A+aa*bb.*cc;end;HH=[10^-5:3*10^-4:10^-3,.0015,.002:.001:.04];figure(1)

36、;contour(r,z,A,HH);xlabel('radical \rho'); ylabel('high z'); colorbar;,5.2.7 柱體內(nèi)溫度場(chǎng)分布之三( I0 的應(yīng)用),半徑為ρ 0高度為L(zhǎng)的圓柱體，側(cè)面有均勻分布的熱流流入，強(qiáng)度為q0, 上、下底溫度保持u0。求圓柱體內(nèi)穩(wěn)定的溫度分布。定解問題為：其解析解為：其中I0為第一類零階修正（或第一類零階虛宗量）貝塞

37、爾函數(shù)。取u0=3; ρ0 =0.6, L=2.5, q0=1用下列程序畫圖，結(jié)果如后頁所示,%ex412; (p120) % 溫度場(chǎng)分布之三；clear; u0=3; r0=0.6; L=2.5; N=21; q0k=1;[r,z]=meshgrid(-r0:0.01:r0,0:0.01:L);A=u0;for n=1:2:N;b=n*pi/L; a=4*q0k/L/b^2; I1=(besseli(0,b

38、*(r0-0.01))-besseli(0,b*r0))/ (b*0.01); c=besseli(0,b*r)/I1;A=A+a*c.*sin(b*z);end;figure(1);contourf(r, z, A, 20);xlabel('radical /\rho'); ylabel('high / z'); colorbar;,5.2.8 柱體內(nèi)溫度場(chǎng)分布之四( I0 的應(yīng)用),半徑

39、為ρ 0高度為L(zhǎng)的圓柱體，上、下底溫度分別保持u1 、 u2，側(cè)面溫度為:求圓柱體內(nèi)穩(wěn)定的溫度分布。定解問題為：其解析解為：取u1=-5; u2=30; ρ0 =0.6, L=2.5,用下列程序畫圖,結(jié)果如后頁所示,%ex413; (p122) % 溫度場(chǎng)分布之四；clear; u1=-5; u2=30; L=2.5; r0=0.6; N=11; x=-r0:0.01:r0; y=0:0.005:L;[r,z]

40、=meshgrid(x,y);u=u1+(u2-u1)/L*z; for n=1:2:N;a=16*u2/(n*pi)^3; b=n*pi/L;u=u-a*besseli(0,b*r)./besseli(0,b*r0).*sin(b*z);end;figure(1);contour(r,z, u, 40);xlabel('radical / \rho');ylabel('high / z'

41、);colorbar;,5.2.9 電子透鏡(J0、I0 的應(yīng)用),如圖所示，電子光學(xué)透鏡的某一部件由兩個(gè)中空?qǐng)A柱筒組成，其電勢(shì)分別為+v0和-v0，在圓柱中間隙縫的邊緣處電勢(shì)可近似為u=v0sin(πz/2δ),求圓柱體內(nèi)的電勢(shì)分布。圓柱體兩端邊界條件可近似為u|=+-l=+-v0, 圓柱體的半徑為ρ0。取柱坐標(biāo)系，由于電勢(shì)為z的奇函數(shù)，可以取通過縫隙中間的中點(diǎn)為柱坐標(biāo)系的原點(diǎn)，則定解問題為：,令：u=v+w,將它分解為兩個(gè)定解

42、問題：其解析解分別為：取u1=-5; u2=30; ρ0 =0.6, L=2.5,用下列程序畫圖,結(jié)果如后頁所示,%ex414; (p125) % 電子透鏡；clear; a=[]; N=50; j=1; v0=100; yy=inline('besselj(0,xx)','xx');for k=1:N; while(besselj(0,j)*besselj(0,j+1)

43、>0), j=j+1;end;%粗估零點(diǎn)；q=fzero(yy,j);a=[a,q];j=j+1; %精估貝塞爾函數(shù)零點(diǎn)；end; h=1.5; r0=0.5; del=0.05; A=0; B=0;[r,z]=meshgrid(0:r0/20:r0,-h:h/20:h);for n=1:N;bb=besselj(0,a(n)*r/r0)/(a(n)*besselj(1,a(n)));cc=sinh(a(n)

44、*z/r0)/sinh(a(n)*h/r0);A=A+bb.*cc; K=n*pi/h;dd=besseli(0,K*r)/besseli(0,K*r0); SS=h^2-(2*n*del)^2;ee=1/(n*pi)*((-1)^(n+1)+h^2/(SS+eps)*cos(K*del))*sin(K*z);B=B+dd.*ee;end; C=2*v0*(A+B); m=min(min(C)); M=max(max(

45、C));figure(1);surfc(r,z,C);axis([0,r0,-h,h,1.1*m,1.1*M]);,5.2.10 柱體外的電勢(shì)( K0 的應(yīng)用),半徑為ρ 0高度為L(zhǎng)的導(dǎo)體圓柱殼，上下底和側(cè)面被絕緣物質(zhì)所隔離，柱殼側(cè)面的電勢(shì)為z，上下底接地。求柱體外的電勢(shì)分布。定解問題為：其解析解為：取ρ 0 =0.5, L=3,用下列程序畫圖,結(jié)果如后頁所示,%ex415; (p128) % 柱體外電勢(shì)；cle

46、ar; a=[]; N=50; j=1; v0=100; h=3; r0=0.5; del=0.05;[r,z]=meshgrid(r0:r0/10:6*r0,0:h/60:h);A=0; for n=1:N; K=n*pi/h;bb=besselk(0,K*r)/besselk(1,K*r0).*sin(K*z);A=A+(-1)^n*2/K*bb; end;figure(1);contour(z,r

47、,A,20);axis([0,h,r0,5*r0]);ylabel('radical /r'); xlabel('high /z'); colorbar;,5.3 泊松方程與格林函數(shù),5.3.1 矩域的泊松方程,在矩形區(qū)域上 求解: 且在邊界上的值為零，定界問題為：,其解為

48、：,%ex416; (p130) % 泊松方程之解clear; a=5; b=5; N=20;x=0:0.1:a; y=-b/2:0.1:b/2; [X,Y]=meshgrid(x,y);Z=X.*Y.*(a^3-X.^3)/12;for n=1:N; K=n*pi/a;aa=(-1)^n*(n*pi)^2+2*(1-(-1)^n);bb=b/a*sin(K*X).*sinh(K*Y)/(K^5*sinh(K*b/

49、2));Z=Z+aa*bb;end;figure(1); mesh(X,Y,Z); colormap(hot); view(119,7);,上圖為解析解取前20項(xiàng)所得的結(jié)果。下圖是采用pdetool所得的結(jié)果。邊界條件取為齊次Dirichlet源項(xiàng)取為f=x.^2.*y,5.3.2 球域的泊松方程與格林函數(shù),在半徑為a的導(dǎo)體球內(nèi)(或球外)，距球心為r0處放置電荷量為4πε0q的點(diǎn)電荷, 求它形成的靜電場(chǎng)。

50、定界問題為：,由鏡像法不難求得其解為：,%ex417; (p133) % 球域格林函數(shù)；clear; a=1; r0=2; r1=a^2/r0; q=1; q1=-q*a/r0; [x,y]=meshgrid(-3:0.01:3,-1:.01:4);[theta,r]=cart2pol(y,x); r(r<=1)=NaN; r(r==2)=NaN;u1=q./sqrt(r0^2-2*r0*r.*cos

51、(theta)+r.^2);u2=q1./sqrt(r1^2-2*r1*r.*cos(theta)+r.^2);u=u1+u2; t=pi*(0:1/20:2); sx=0.1*cos(t); sy=r0+0.1*sin(t);figure(1);contour(x,y,u,[-1:20,3000]);hold on;box on;colorbar;[Ex,Ey]=gradient(-u); streamline(x,y,E

52、x,Ey,sx,sy); axis equal;plot(a*exp(i*t),'r');title('球域外點(diǎn)電荷的電場(chǎng)');figure(2);contour(x,y,u1,[-1:20,3000]);hold on;box on;colorbar;[Ex,Ey]=gradient(-u1); streamline(x,y,Ex,Ey,sx,sy); axis equal;plot(a*exp

53、(i*t),'r');title('球域外源電荷的電場(chǎng)');figure(3);contour(x,y,u2,11);hold on;box on;colorbar;[Ex,Ey]=gradient(-u2); streamline(x,y,Ex,Ey,sx,sy); axis equal;plot(a*exp(i*t),'r');title('球域外像電荷的電場(chǎng)')

54、;,%ex417; (p133) % 球域格林函數(shù)；clear; a=2; r0=1; r1=a^2/r0; q=1; q1=-q*a/r0; [x,y]=meshgrid(-a:0.01:a, -a:0.01:a);[theta,r]=cart2pol(y,x); r(r>=a)=NaN; u1=q./(sqrt(r0^2-2*r0*r.*cos(theta)+r.^2)+eps);u2=q1./sqrt(r1

55、^2-2*r1*r.*cos(theta)+r.^2);u=u1+u2; t=pi*(0:1/20:2); sx=0.1*cos(t); sy=r0+0.1*sin(t);figure(1);contour(x,y,u,[-1:20,3000]);hold on;box on;colorbar;[Ex,Ey]=gradient(-u);streamline(x,y,Ex,Ey,sx,sy);axis equal;plot(a

56、*exp(i*t),‘r’);axis([-a,a,-a,a]);title(‘球域內(nèi)點(diǎn)電荷的電場(chǎng)');figure(2);contour(x,y,u1,[-1:20,3000]);hold on;box on;colorbar;[Ex,Ey]=gradient(-u1);streamline(x,y,Ex,Ey,sx,sy);axis equal;plot(a*exp(i*t),‘r’);axis([-a,a,-a,a]

57、);title(‘球域內(nèi)源電荷的電場(chǎng)');figure(3);contour(x,y,u2,11);hold on;box on;colorbar;[Ex,Ey]=gradient(-u2);streamline(x,y,Ex,Ey,sx,sy);axis equal;plot(a*exp(i*t),‘r’);axis([-a,a,-a,a]);title(‘球域內(nèi)像電荷的電場(chǎng)');,5.3.3 圓域的泊松方程與格

58、林函數(shù),在半徑為a的導(dǎo)體圓柱外，距軸心為r0處放置單位長(zhǎng)度電量為λ =ε0q的線電荷, 求它形成的靜電場(chǎng)。 定界問題為：,其解為：其中：r0為源點(diǎn)位置，r為場(chǎng)點(diǎn)位置,%ex419; (p136) % 圓域格林函數(shù)；clear; a=1; r0=2; r1=a^2/r0; q=1; [x,y]=meshgrid(-5:0.02:5,-5:.02:5); [theta,r]=cart2pol(y,x); r(

59、r<=a)=NaN; r(r==2)=NaN;u1=q/(2*pi)*log(a./sqrt(r0^2-2*r0*r.*cos(theta)+r.^2));u2=-q/(2*pi)*log(r1./sqrt(r1^2-2*r1*r.*cos(theta)+r.^2));u=u1+u2; t=pi*(0:1/20:2); sx=0.1*cos(t);sy=r0+0.1*sin(t);figure(1);contour(x,y

60、,u,15);hold on;box on;colorbar; [Ex,Ey]=gradient(-u); streamline(x,y,Ex,Ey,sx,sy);axis equal; plot(a*exp(i*t),'r'); axis([-5,5,-5,5]); title('圓域外線電荷的電場(chǎng)');figure(2);contour(x,y,u1,15);hold on;box on;c

61、olorbar; [Ex,Ey]=gradient(-u1);streamline(x,y,Ex,Ey,sx,sy);axis equal; plot(a*exp(i*t),'r');axis([-5,5,-5,5]);title('圓域外源電荷的電場(chǎng)');figure(3);contour(x,y,u2,7);hold on;box on;colorbar; [Ex,Ey]=gradien

62、t(-u2);streamline(x,y,Ex,Ey,sx,sy);axis equal; plot(a*exp(i*t),'r');axis([-5,5,-5,5]);title('圓域外像電荷的電場(chǎng)');,%ex420; (p137) % 圓域格林函數(shù)；clear; a=2; r0=1; r1=a^2/r0; q=1; [x,y]=meshgrid(-a:0.01:a,-a:.01:a

63、);[theta,r]=cart2pol(y,x); r(r>=a)=NaN; r(r==r0)=NaN; u1=q/(2*pi)*log(a./sqrt(r0^2-2*r0*r.*cos(theta)+r.^2));u2=-q/(2*pi)*log(r1./sqrt(r1^2-2*r1*r.*cos(theta)+r.^2));u=u1+u2;t1=pi*(0:1/10:2);t=pi*(0:1/10:2);sx=0.

64、1*cos(t);sy=r0+0.1*sin(t);figure(1);contour(x,y,u,15);hold on;box on;colorbar; [Ex,Ey]=gradient(-u);streamline(x,y,Ex,Ey,sx,sy);axis equal; plot(a*exp(i*t1),'r');axis([-a,a,-a,a]);title('圓域內(nèi)線電荷的電場(chǎng)');f

65、igure(2);contour(x,y,u1,15);hold on;box on;colorbar; [Ex,Ey]=gradient(-u1);streamline(x,y,Ex,Ey,sx,sy);axis equal; plot(a*exp(i*t1),'r');axis([-a,a,-a,a]);title('圓域內(nèi)源電荷的電場(chǎng)');figure(3);contour(x,y,u2,7)