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

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

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

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

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

5、代入(20)將A和B代入（9）有、Matlab程序?qū)崿F(xiàn)2.1 當(dāng)右邊界條件為時 用Matlab對方程在邊界條件下根據(jù)a的變化來分析方程的解程序代碼a=-10,-6,-4,0,6,10,u=40,R=13; for i=1:6 title('Bessel函數(shù)在a變化時的情形') 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變化時y隨x的增加具有衰減振蕩性，即y的振幅隨x值增大而減小

8、。2、 并且方程解y隨x的變化是成周期性的，變化的周期T的大小跟a的大小無關(guān)。3、 方程（1）在a變化時，方程解y隨x的變化圖像有固定交點，并且交點位置與a的大小無關(guān)。用Matlab對方程在邊界條件下根據(jù)R的變化來分析方程的解程序代碼a=6,u=50,R=30,40,100,130,170,200; for i=1:6 title('bessel函數(shù)在R變化時的情形') 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變化時y隨x的增加具有衰減振蕩性，即y的振幅隨x值增大而減小。2、并且方程解y隨x變化是成周期性的，變化的周期T的大小跟R的大小無關(guān)。3、方程（1）在R變化時，方程解y隨x的變化圖像有固定交點，并且交點位置與R的取值無關(guān)。用Matlab對方程在邊界條件下根據(jù)的變化來分析方程的解程序代碼a=33,R=135,u=30,40,50,70,80,100; for i=1:6 title('bessel函數(shù)在u變化時的情形') 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）在變化時y隨x的增加具有衰減振蕩性，即y的振幅隨x值增大而減小。2、 并且方程解y隨x的變化是成周期性，變化的周期T的大小跟的大小無關(guān)。3、 方程（1）在變化時，方程解y隨x的變化圖像無固定交點，可知交點位置與值有關(guān)。2.2 當(dāng)右邊界條件是Y(R)=0時 用Matlab對方程在邊界條件下根據(jù)a的變化來分析方程的解程序代碼a=-10,-6,-4,0,6,10,u=40,R=13; for i=1:6 title('bessel函數(shù)在a變化時的情形') 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變化時y隨x的增加具有衰減振蕩性，即y的振幅隨x值增大而減小。2、 并且方程解y隨x的變化是成周期性，變化的周期T的大小跟a的大小無關(guān)。3、 方程（1）在a變化時，方程解y隨x的變化圖像有固定交點，并且交點位置與a的大小無關(guān)。2.2.2 用Matlab對方程在邊界條件下根據(jù)R的變化來分析方程的解程序代碼a=6,u=50,R=30,40,100,130,170,200; for i=1:

18、6 title('bessel函數(shù)在R變化時的情形') 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變化時y隨x的增加具有衰減振蕩性，即y的振幅隨x值增大而減小。2、 并且方程解y隨x的變化是成周期性的，變化的周期T的大小跟R的大小無關(guān)。3、 方程（1）在R變化時，方程解y隨x的變化圖像有固定交點，可知交點位置與R值無關(guān)。2.2.3 用Matlab對方程在邊界條件下根據(jù)的變化來

21、分析方程的解程序代碼a=33,R=135,u=30,40,50,70,80,100; for i=1:6 title('bessel函數(shù)在u變化時的情形') 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',