的特解進(jìn)行討論并用Matlab程序?qū)崿F(xiàn)繪制特解的變化圖像

1、、前言： 到目前為止已經(jīng)有不少人對(duì)Bessel方程進(jìn)行了討論，但對(duì)Bessel方程的特定邊界條件下的特解情況的討論卻很少。本文主要是針對(duì)一類零階Bessel方程在如下邊界條件的特解進(jìn)行討論并用Matlab程序?qū)崿F(xiàn)繪制特解的變化圖像。Matlab是類似于Fortran和C的一種語(yǔ)言，雖然很難簡(jiǎn)短描述其特點(diǎn)，但對(duì)于科學(xué)計(jì)算而言，其主要特色有：1) 易于編程2) 整數(shù)、實(shí)數(shù)、復(fù)數(shù)之間的統(tǒng)一性3) 高精度及擴(kuò)充的數(shù)值范圍4) 綜合性的數(shù)學(xué)工具庫(kù)5) 包含圖形用戶界面在內(nèi)的功能完備的圖形工具6) 與傳統(tǒng)編程語(yǔ)言的接口7) Matlab程序的可移植性本文便是利用的強(qiáng)大圖形工具繪制出方程特解的圖像，以便形象

2、地描述Bessel方程特解的振蕩性、衰減性、周期性等。、數(shù)學(xué)模型的求解及分析： 數(shù)學(xué)模型：1.1 方程通解的求解過(guò)程 首先對(duì)方程 (1)進(jìn)行化簡(jiǎn)。引入，則有 = (2)=將方程(1)兩邊同時(shí)乘以后，經(jīng)整理后得: (3)將方程(2)代入到方程(3)有: (4)對(duì)比Bessel標(biāo)準(zhǔn)方程（這里是非負(fù)常數(shù)，不一定是正整數(shù)）知，方程（4）其實(shí)就是一個(gè)的Bessel方程；即零階的Bessel方程。根據(jù)參考文獻(xiàn)158頁(yè)定理11:若方程 中系數(shù)具有這樣的性質(zhì),即和均能展成的冪級(jí)數(shù),且收斂區(qū)間為,則方程有形如 即 的特解,這里,是一個(gè)特定的常數(shù)。級(jí)數(shù)也以為收斂區(qū)間。將方程（4）改寫(xiě)成 易見(jiàn)，它滿足定理11的條件

3、，且, , 按展開(kāi)成冪級(jí)數(shù),它的收斂區(qū)間為，由定理11，方程有形如 (5)的解.這里。而和是待定常數(shù)。將（5）代如（4）有 按t的同冪次項(xiàng)整理有：令各項(xiàng)系數(shù)等于零，得一系列的代數(shù)方程： (6)因?yàn)?，故?6)的第一個(gè)方程解得：考慮時(shí)方程(4)的一個(gè)特解.這時(shí)總可以從(6)中逐個(gè)地確定所有的系數(shù)。把代入（6）得到 k=2,3,對(duì)討論奇偶兩種情況，分別有 k=1,2,3,. 從而求得 k=1,2,3,. k=1,2,3,將各代入方程(5)得到方程(4)的一個(gè)解 (7)不妨令 <2>則(7)變?yōu)?根據(jù)參考文獻(xiàn)P404 (7.9): (n=0時(shí),右邊第二項(xiàng)不存在)可以得到 (8)當(dāng)為實(shí)數(shù)時(shí)

4、, 是一個(gè)衰減振蕩函數(shù), 在有奇性。 、都是Bessel方程的解，且其中任意兩個(gè)都是線性無(wú)關(guān)的。所以Bessel方程（4）的通解可寫(xiě)成： (9)1.2 在左右邊界條件下方程特解的求解過(guò)程1.2.1左邊界條件根據(jù)Bessel函數(shù)性質(zhì)有 即有 (12)現(xiàn)在將(8)和(12)代入(11)有 (13) 當(dāng)右邊界條件為時(shí)有 =0 (14)將(13)和(14)組成方程組聯(lián)合求解和有 (15)代入（13）得： 即將(16)代入(15)有 (17)將A和B代入（9）有當(dāng)右邊界條件是Y(R)=0時(shí)有令有 并且 (18)將方程（13）和（18）組成方程組求解A和B由（18）得： (19)代入(13)有即將(19)

5、代入(20)將A和B代入（9）有、Matlab程序?qū)崿F(xiàn)2.1 當(dāng)右邊界條件為時(shí) 用Matlab對(duì)方程在邊界條件下根據(jù)a的變化來(lái)分析方程的解程序代碼a=-10,-6,-4,0,6,10,u=40,R=13; for i=1:6 title('Bessel函數(shù)在a變化時(shí)的情形') xlabel('自變量X') ylabel('函數(shù)值Y') b=1-a(i);t=sqrt(u); Q=60;B=besselj(1,t*R)*Q/(besselj(1,t*R)*(a(i)*bessely(0,t)-b*t*bessely(1,t)-bessely(1,t

6、*R)*(a(i)*besselj(0,t)-b*t*besselj(1,t);A=-bessely(1,t*R)*Q/(besselj(1,t*R)*(a(i)*bessely(0,t)-b*t*bessely(1,t)-bessely(1,t*R)*(a(i)*besselj(0,t)-b*t*besselj(1,t);x=0:0.05:50;y=A*besselj(0,x)+B*bessely(0,x);switch icase 1,plot(x,y,'r-');case 2,plot(x,y,'g:');case 3,plot(x,y,'b-.&

7、#39;); case 4,plot(x,y,'y-');case 5,plot(x,y,'m:');case 6,plot(x,y,'c-.'); end;hold on;end;hold off;legend ('a=-10','a=-6','a=-4','a=0','a=6','a=10',0);形成方程解的變化圖像:圖像分析1、 從方程解的變化圖像中很容易看出Bessel方程（1）在a變化時(shí)y隨x的增加具有衰減振蕩性，即y的振幅隨x值增大而減小

8、。2、 并且方程解y隨x的變化是成周期性的，變化的周期T的大小跟a的大小無(wú)關(guān)。3、 方程（1）在a變化時(shí)，方程解y隨x的變化圖像有固定交點(diǎn)，并且交點(diǎn)位置與a的大小無(wú)關(guān)。用Matlab對(duì)方程在邊界條件下根據(jù)R的變化來(lái)分析方程的解程序代碼a=6,u=50,R=30,40,100,130,170,200; for i=1:6 title('bessel函數(shù)在R變化時(shí)的情形') xlabel('自變量X') ylabel('函數(shù)值Y') b=1-a;t=sqrt(u); Q=60;B=besselj(1,t*R(i)*Q/(besselj(1,t*R(i

9、)*(a*bessely(0,t)-b*t*bessely(1,t)-bessely(1,t*R(i)*(a*besselj(0,t)-b*t*besselj(1,t);A=-bessely(1,t*R(i)*Q/(besselj(1,t*R(i)*(a*bessely(0,t)-b*t*bessely(1,t)-bessely(1,t*R(i)*(a*besselj(0,t)-b*t*besselj(1,t);x=0:0.05:50;y=A*besselj(0,x)+B*bessely(0,x);switch icase 1,plot(x,y,'r-');case 2,plo

10、t(x,y,'g:');case 3,plot(x,y,'b-.'); case 4,plot(x,y,'y-');case 5,plot(x,y,'m:');case 6,plot(x,y,'c-.'); end;hold on;end;hold off;legend ('R=30','R=40','R=100','R=130','R=170','R=200',0)形成方程解的變化圖像:圖像分析1、從方程解的變化圖像中

11、很容易看出Bessel方程（1）在R變化時(shí)y隨x的增加具有衰減振蕩性，即y的振幅隨x值增大而減小。2、并且方程解y隨x變化是成周期性的，變化的周期T的大小跟R的大小無(wú)關(guān)。3、方程（1）在R變化時(shí)，方程解y隨x的變化圖像有固定交點(diǎn)，并且交點(diǎn)位置與R的取值無(wú)關(guān)。用Matlab對(duì)方程在邊界條件下根據(jù)的變化來(lái)分析方程的解程序代碼a=33,R=135,u=30,40,50,70,80,100; for i=1:6 title('bessel函數(shù)在u變化時(shí)的情形') xlabel('自變量X') ylabel('函數(shù)值Y') b=1-a; t=sqrt(u(

12、i); Q=60;B=besselj(1,t*R)*Q/(besselj(1,t*R)*(a*bessely(0,t)-b*t*bessely(1,t)-bessely(1,t*R)*(a*besselj(0,t)-b*t*besselj(1,t);A=-bessely(1,t*R)*Q/(besselj(1,t*R)*(a*bessely(0,t)-b*t*bessely(1,t)-bessely(1,t*R)*(a*besselj(0,t)-b*t*besselj(1,t);x=0:0.05:50;y=A*besselj(0,x)+B*bessely(0,x);switch icase 1

13、,plot(x,y,'r-');case 2,plot(x,y,'g:');case 3,plot(x,y,'b-.'); case 4,plot(x,y,'y-');case 5,plot(x,y,'m:');case 6,plot(x,y,'c-.'); end;hold on;end;hold off;legend ('u=30','u=40','u=50','u=70','u=80','u=100'

14、;,0)形成方程解的變化圖像:圖像分析1、 從方程解的變化圖像中很容易看出Bessel方程（1）在變化時(shí)y隨x的增加具有衰減振蕩性，即y的振幅隨x值增大而減小。2、 并且方程解y隨x的變化是成周期性，變化的周期T的大小跟的大小無(wú)關(guān)。3、 方程（1）在變化時(shí)，方程解y隨x的變化圖像無(wú)固定交點(diǎn)，可知交點(diǎn)位置與值有關(guān)。2.2 當(dāng)右邊界條件是Y(R)=0時(shí) 用Matlab對(duì)方程在邊界條件下根據(jù)a的變化來(lái)分析方程的解程序代碼a=-10,-6,-4,0,6,10,u=40,R=13; for i=1:6 title('bessel函數(shù)在a變化時(shí)的情形') xlabel('自變量X&

15、#39;) ylabel('函數(shù)值Y') b=1-a(i);t=sqrt(u); Q=60;B=besselj(0,t*R)*Q/(besselj(0,t*R)*(a(i)*bessely(0,t)-b*t*bessely(1,t)-bessely(0,t*R)*(a(i)*besselj(0,t)-b*t*besselj(1,t);A=-bessely(0,t*R)*Q/(besselj(0,t*R)*(a(i)*bessely(0,t)-b*t*bessely(1,t)-bessely(0,t*R)*(a(i)*besselj(0,t)-b*t*besselj(1,t);x

16、=0:0.05:50;y=A*besselj(0,x)+B*bessely(0,x);switch icase 1,plot(x,y,'r-');case 2,plot(x,y,'g:');case 3,plot(x,y,'b-.'); case 4,plot(x,y,'y-');case 5,plot(x,y,'m:');case 6,plot(x,y,'c-.'); end;hold on;end;hold off;legend ('a=-10','a=-6',&

17、#39;a=-4','a=0','a=6','a=10',0);形成方程解的變化圖像:分析圖像1、 從方程解的變化圖像中很容易看出Bessel方程（1）在a變化時(shí)y隨x的增加具有衰減振蕩性，即y的振幅隨x值增大而減小。2、 并且方程解y隨x的變化是成周期性，變化的周期T的大小跟a的大小無(wú)關(guān)。3、 方程（1）在a變化時(shí)，方程解y隨x的變化圖像有固定交點(diǎn)，并且交點(diǎn)位置與a的大小無(wú)關(guān)。2.2.2 用Matlab對(duì)方程在邊界條件下根據(jù)R的變化來(lái)分析方程的解程序代碼a=6,u=50,R=30,40,100,130,170,200; for i=1:

18、6 title('bessel函數(shù)在R變化時(shí)的情形') xlabel('自變量X') ylabel('函數(shù)值Y') b=1-a;t=sqrt(u); Q=60;B=besselj(0,t*R(i)*Q/(besselj(0,t*R(i)*(a*bessely(0,t)-b*t*bessely(1,t)-bessely(0,t*R(i)*(a*besselj(0,t)-b*t*besselj(1,t);A=-bessely(0,t*R(i)*Q/(besselj(0,t*R(i)*(a*bessely(0,t)-b*t*bessely(1,t)-b

19、essely(0,t*R(i)*(a*besselj(0,t)-b*t*besselj(1,t);x=0:0.05:50;y=A*besselj(0,x)+B*bessely(0,x);switch icase 1,plot(x,y,'r-');case 2,plot(x,y,'g:');case 3,plot(x,y,'b-.'); case 4,plot(x,y,'y-');case 5,plot(x,y,'m:');case 6,plot(x,y,'c-.'); end;hold on;end

20、;hold off;legend ('R=30','R=40','R=100','R=130','R=170','R=200',0)形成方程解的變化圖像:分析圖像1、 從方程解的變化圖像中很容易看出Bessel方程（1）在R變化時(shí)y隨x的增加具有衰減振蕩性，即y的振幅隨x值增大而減小。2、 并且方程解y隨x的變化是成周期性的，變化的周期T的大小跟R的大小無(wú)關(guān)。3、 方程（1）在R變化時(shí)，方程解y隨x的變化圖像有固定交點(diǎn)，可知交點(diǎn)位置與R值無(wú)關(guān)。2.2.3 用Matlab對(duì)方程在邊界條件下根據(jù)的變化來(lái)

21、分析方程的解程序代碼a=33,R=135,u=30,40,50,70,80,100; for i=1:6 title('bessel函數(shù)在u變化時(shí)的情形') xlabel('自變量X') ylabel('函數(shù)值Y') b=1-a; t=sqrt(u(i); Q=60;B=besselj(0,t*R)*Q/(besselj(0,t*R)*(a*bessely(0,t)-b*t*bessely(1,t)-bessely(0,t*R)*(a*besselj(0,t)-b*t*besselj(1,t);A=-bessely(0,t*R)*Q/(besse

22、lj(0,t*R)*(a*bessely(0,t)-b*t*bessely(1,t)-bessely(0,t*R)*(a*besselj(0,t)-b*t*besselj(1,t);x=0:0.05:50;y=A*besselj(0,x)+B*bessely(0,x);switch icase 1,plot(x,y,'r-');case 2,plot(x,y,'g:');case 3,plot(x,y,'b-.'); case 4,plot(x,y,'y-');case 5,plot(x,y,'m:');case 6,plot(x,y,'c-.'); end;hold on;end;hold off;legend ('u=30','u=40','u=50',