龍姝明， 冉啟武， 董蓓蓓

(1.陜西理工學(xué)院 物理與電信工程學(xué)院， 陜西 漢中 723000；2.陜西理工學(xué)院 電氣工程學(xué)院， 陜西 漢中 723003)

0 引 言

電器產(chǎn)品越人性化，電器設(shè)備的電路就越復(fù)雜。要計(jì)算出滿足用戶需求的電器設(shè)備電路元件參數(shù)，必須求出電路所有支路電流電壓的演化規(guī)律，手工幾乎是不可能快速解決這一問(wèn)題的。編程列方程和解方程是行之有效的方法。

對(duì)于復(fù)雜多輸入多輸出動(dòng)態(tài)電路，要得到全部支路電流電壓演化規(guī)律，面臨兩個(gè)難題。一是列寫(xiě)時(shí)域方程比較困難，因?yàn)樵碴P(guān)系不再是簡(jiǎn)單代數(shù)關(guān)系，支路電壓與支路電流的關(guān)系可能是微分積分關(guān)系，不再是簡(jiǎn)單的一個(gè)輸入一個(gè)輸出的電路系統(tǒng)。因而列出的時(shí)域方程只能是微分積分方程組[1]，難以簡(jiǎn)化成一個(gè)高階微分方程或一階狀態(tài)方程組。二是列出的微分積分方程組，在求解過(guò)程中必須面對(duì)手工無(wú)法解決的高次代數(shù)方程求根問(wèn)題。解決上述第一個(gè)難題的有效方法是，用Laplace積分變換將時(shí)域動(dòng)態(tài)電路問(wèn)題轉(zhuǎn)化為復(fù)頻域(即Laplace變換的像域，又稱為s域)的“像電路”問(wèn)題[2]，再設(shè)計(jì)程序自動(dòng)列寫(xiě)像電路節(jié)點(diǎn)電壓方程或回路電流方程并解方程。解決上述第二個(gè)難題的唯一方法是，用軟件(Mathematica，Maple，Matlab等)編程來(lái)分解高次代數(shù)多項(xiàng)式為多個(gè)簡(jiǎn)單因子的乘積，為求Laplace逆積分變換給出所有支路電流電壓時(shí)域演化規(guī)律奠定基礎(chǔ)。

自動(dòng)列寫(xiě)復(fù)雜多輸入多輸出動(dòng)態(tài)電路方程，快捷求解全部支路電流電壓隨時(shí)間演化的規(guī)律，直觀表示各支路電流電壓演化特性，用能量守恒來(lái)檢驗(yàn)求解結(jié)果的正確性，用程序來(lái)實(shí)現(xiàn)這一目標(biāo)，對(duì)于復(fù)雜動(dòng)態(tài)電路而言，是一項(xiàng)有挑戰(zhàn)性的工作。

1 拉普拉斯積分變換在動(dòng)態(tài)電路求解中的作用

在我們研究的多輸入多輸出動(dòng)態(tài)電路中，先假設(shè)：①時(shí)域信號(hào)源函數(shù)總可以用一組帶單位階躍函數(shù)因子的數(shù)學(xué)函數(shù)表示，并且這些函數(shù)在Laplace積分變換下的像函數(shù)很容易手工或編程求得；②電路所有動(dòng)態(tài)元件積累的初始能量為零(只關(guān)注長(zhǎng)時(shí)間后電路的穩(wěn)態(tài)響應(yīng))；③電路中所有信號(hào)源都是內(nèi)阻為零的理想電壓源或內(nèi)阻為無(wú)限大的理想電流源(如果內(nèi)阻有限非零且已知，則可以將其歸并到相應(yīng)電阻元件中)。

取電感元件初始電流和電容元件初始電壓為0，將電阻R、電容C、電感L時(shí)域電流電壓約束關(guān)系(即伏安關(guān)系)取Laplace變換映射到像(s)域，則有表1第四列的s域伏安關(guān)系。

表1 電阻、電容、電感時(shí)域、像域伏安關(guān)系即像域電阻模型

比較元件時(shí)域與像域伏安關(guān)系發(fā)現(xiàn)，如果將電感和電容看成導(dǎo)納分別為GL=1/(sL)和GC=sC的電導(dǎo)元件，那么時(shí)域動(dòng)態(tài)電路在Laplace變換的像域(s域)中，可看成電阻電路(參見(jiàn)表1)。于是不再需要列寫(xiě)動(dòng)態(tài)電路的時(shí)域微積分方程，只要列出各擴(kuò)展支路像電壓像電流滿足的代數(shù)方程[3]，就可以解出全部支路的像電壓像電流的s域函數(shù)。再取逆Laplace積分變換(也即依據(jù)黎曼梅林公式求留數(shù)和)就能給出全部支路的時(shí)域電壓電流函數(shù)。

結(jié)論是，只要取Laplace積分變換，解復(fù)雜動(dòng)態(tài)電路的問(wèn)題就不再是難題。

2 s域擴(kuò)展支路模型

節(jié)點(diǎn)電壓分析方法中，不允許出現(xiàn)理想電壓源支路。要允許電路中既有理想電流源又有理想電壓 源存在，需要擴(kuò)展像域支路的概念。我們?cè)诶硐腚妷涸粗分写搿疤摂M小電阻r”(求出s域節(jié)點(diǎn)像電壓后再取r趨近于零的極限，就可以消去r)，將一條單一元件支路擴(kuò)展為三個(gè)元件組成的一條擴(kuò)展支路，以便我們的求解程序可以處理有理想電壓源支路的情況。我們的一條擴(kuò)展支路是三個(gè)元件的有序組合，即一個(gè)理想電壓源Us(s)串聯(lián)一個(gè)導(dǎo)納元件G(s)，再并聯(lián)一個(gè)電流源Is(s)，擴(kuò)展支路如圖1，其中導(dǎo)納元件可以是電阻、電容、電感，也可以是這三種元件的串、并聯(lián)組合。

圖1 復(fù)雜動(dòng)態(tài)電路分析中 的s域擴(kuò)展支路模型

設(shè)Us(s)=LT[us(t)]，Is(s)=LT[is(t)]分別表示對(duì)us(t)，is(t)取Laplace積分變換。約定實(shí)際支路與擴(kuò)展支路的對(duì)應(yīng)關(guān)系用表2描述。引入擴(kuò)展支路模型，為編程自動(dòng)列寫(xiě)復(fù)雜動(dòng)態(tài)電路s域節(jié)點(diǎn)電壓方程組做好了準(zhǔn)備。

表2 實(shí)際元件支路與擴(kuò)展支路模型參數(shù)對(duì)應(yīng)關(guān)系

3 建立描述電路拓?fù)潢P(guān)系的重要矩陣

為了最大限度簡(jiǎn)化建立電路方程的過(guò)程，首要問(wèn)題是電路拓?fù)潢P(guān)系和元件參數(shù)的輸入與表示。這需要將電路擴(kuò)展支路和獨(dú)立節(jié)點(diǎn)有序化，即分別編號(hào)索引。方便的做法是,給每個(gè)擴(kuò)展支路設(shè)一個(gè)電流流向,選支路電流流入端為支路電壓的“+”端，流出端為支路電壓的“-”端。這樣選擇的結(jié)果是電路中所有元件消耗的功率總和為零(即能量守恒)，信號(hào)源消耗功率小于0時(shí)，表示信號(hào)源為電路提供大于0的功率。

3.1 輸入時(shí)域電壓源矢量usv[t]和電流源矢量isv[t]

約定支路電流流入電壓源的“+”端時(shí)，電壓源電壓取正，否則，電壓源電壓取負(fù)。約定電流源電流與支路電流流向一致時(shí)，電流源電流取正，否則，電流源電流取負(fù)。由電流源和電壓源函數(shù)乘正負(fù)號(hào)來(lái)構(gòu)成時(shí)域電壓源矢量usv[t]和電流源矢量isv[t]。再取Laplace變換，給出s域電壓源矢量Usv和電流源矢量Isv。

3.2 輸入電路元件參數(shù)及拓?fù)潢P(guān)系表示矩陣EPTR

設(shè)我們所要求解的電路有n個(gè)節(jié)點(diǎn)，k條擴(kuò)展支路。要最大限度降低表示、輸入電路元件參數(shù)、拓?fù)潢P(guān)系的難度，可以先將電路的k條擴(kuò)展支路編號(hào)為1到k，將n個(gè)節(jié)點(diǎn)編號(hào)為1到n，連接支路最多的節(jié)點(diǎn)賦以最大節(jié)點(diǎn)編號(hào)n。以擴(kuò)展支路為單元輸入電路元件參數(shù)、信號(hào)源像電流Is(s)、像電壓Us(s)及電路拓?fù)潢P(guān)系信息。設(shè)所要研究求解的電路元件參數(shù)及拓?fù)潢P(guān)系表示矩陣為EPTR(Element Parameters and Topology Relations)。規(guī)定EPTR矩陣有k行6列，第j行的6個(gè)元素{bj，intoj，outj，Gj，Usj，Isj}描述第j條擴(kuò)展支路的元件參數(shù)、電壓源和電流源s域像函數(shù)及其與其它支路的拓?fù)潢P(guān)系。b表示擴(kuò)展支路編號(hào)，into和out分別表示擴(kuò)展支路電流流入端節(jié)點(diǎn)和流出端節(jié)點(diǎn)的編號(hào)，G表示擴(kuò)展支路的等效像導(dǎo)納，Us表示擴(kuò)展支路電壓源s域像電壓函數(shù)，Is表示擴(kuò)展支路電流源s域像電流函數(shù)(若Us或Is為負(fù)值，表示該支路中電壓源極性或電流源電流方向與參考方向相反)。EPTR矩陣可表示為

(1)

按照上述規(guī)定對(duì)電路有序化，列出EPTR矩陣后，就可以方便地將電路元件相關(guān)參數(shù)和電路拓?fù)潢P(guān)系清楚地采集并列寫(xiě)出來(lái)，為后續(xù)的計(jì)算機(jī)編程提供方便。

3.3 建立拓?fù)潢P(guān)系矩陣

運(yùn)行程序：

Aa=ConstantArray[0,{n,k}];

For[j=1,j<=k, j++, i=EPTR[[j,2]];Aa[[i, j]]=1;q=EPTR[[j,3]];Aa[[q, j]]=-1];

A=Take[Aa,n-1];

可由EPTR矩陣快速求出拓?fù)潢P(guān)系矩陣A，矩陣A有n-1行(分別表示n-1個(gè)獨(dú)立節(jié)點(diǎn))k列(表示k條支路)。矩陣A描述了每條擴(kuò)展支路關(guān)聯(lián)的電流流入節(jié)點(diǎn)(用1表示)和電流流出節(jié)點(diǎn)(用-1表示)。

3.4 建立電路電導(dǎo)矩陣Gm

用節(jié)點(diǎn)電壓法列寫(xiě)電路方程時(shí)，要知道每個(gè)節(jié)點(diǎn)的自電導(dǎo)和互電導(dǎo)，因此建立節(jié)點(diǎn)電壓方程時(shí)必須列寫(xiě)n-1行n-1列的電導(dǎo)矩陣Gm，它的主對(duì)角線元素是各獨(dú)立節(jié)點(diǎn)的自電導(dǎo)，非主對(duì)角元素是獨(dú)立節(jié)點(diǎn)間的互電導(dǎo)，互電導(dǎo)帶負(fù)號(hào)。我們先建立b行b列的對(duì)角電導(dǎo)矩陣Gdm，它的主對(duì)角線元素對(duì)應(yīng)每條支路的等效電導(dǎo)，非對(duì)角元素均為零。再用矩陣公式

AT=Transpose[A]，Gm=A·Gdm·AT，

(2)

求出Gm矩陣，其中AT是A的轉(zhuǎn)置矩陣。從EPTR中自動(dòng)求出Gm矩陣的程序?yàn)椋?/p>

Gd=ConstantArray[0,{n,k}]；

For[j=1, j<=k, j++,Gdm[[j, j]]=EPTR[[j，4]]]; Gm=A·Gdm·AT。

3.5 節(jié)點(diǎn)的s域電流源矢量Isv

在s域中，流入每一支路的電流源可以集合為一個(gè)電流源矢量Isv，每一支路的電壓源可以集合為一個(gè)電壓源矢量Usv。運(yùn)行程序：

Isv=EPTR[[All,6]];Usv=EPTR[[All,5]];

可以求出Usv及Isv?；蛘撸瑢?duì)usv[t]和isv[t]取Laplace積分變換給出Usv和Isv。再通過(guò)支路電導(dǎo)對(duì)角矩陣Gdm和電路拓?fù)潢P(guān)系矩陣A，可解得流入每個(gè)節(jié)點(diǎn)的總像電流源矢量：

Insv=A·(Gdm·Usv-Isv)。

(3)

3.6 僅有理想電壓源元件的支路索引矩陣ivsm

必須記錄單獨(dú)的理想電壓源構(gòu)成的支路，為求單獨(dú)理想電壓源支路的電流和電壓提供必要的信息。記錄格式為：

ivsm={{單獨(dú)理想電壓源支路編號(hào)，電流流入端節(jié)點(diǎn)號(hào)},…}。

運(yùn)行程序：

ivsb=Position[EPTR,1/r][[All,1]]; Livsm=Length[ivsb]; ivsm={};

For[j=1, j<=Livsm, j++, jb=ivsb[[j]];AppendTo[ivsm,{jb,EPTR[jb,2]]}]]。

可由EPTR自動(dòng)建立ivsm矩陣。

4 復(fù)雜動(dòng)態(tài)電路節(jié)點(diǎn)像電壓方程的列寫(xiě)和求解

復(fù)雜動(dòng)態(tài)電路節(jié)點(diǎn)電壓法是以獨(dú)立節(jié)點(diǎn)像電壓作為未知函數(shù)，將各獨(dú)立節(jié)點(diǎn)像電壓函數(shù)集合為矢量Unv，應(yīng)用基爾霍夫電流、電壓定律和擴(kuò)展支路伏安關(guān)系，可以給出復(fù)雜動(dòng)態(tài)電路關(guān)于n-1個(gè)獨(dú)立節(jié)點(diǎn)的節(jié)點(diǎn)電壓方程組：

Gm·Unv=Insv。

(4)

(5)

利用電路拓?fù)潢P(guān)系矩陣解出支路像電壓矢量Ubv：

Ubv=AT·Unv。

(6)

利用擴(kuò)展支路電路模型，分別解出支路電導(dǎo)元件兩端的像電壓矢量Uzv、支路電導(dǎo)元件中的電流矢量Izv、支路像電流矢量Ibv：

Uzv=Ubv-Usv,Izv=Gdm·Uzv，Ibv=Izv+Isv。

(7)

由單電壓源支路索引矩陣ivsm，依據(jù)節(jié)點(diǎn)電流定律，可以導(dǎo)出單電壓源支路像電流，所需程序?yàn)椋?/p>

For[j=1, j<=Length[ivsm], j++,ib=ivsm[[j,1]];ic=isvm[[j,2]];cin=Drop[Range[1,k],ib];

Ibv[[ib]]=-A[[ic,ib]]*Sum[A[[ic,p]]*Ibv[[p]],{p,cin}]]。

電導(dǎo)元件消耗(包括電容、電感存儲(chǔ))的總功率必定等于所有電源供給的總功率，即有

p(t)=izv(t)uzv(t)+izv(t)usv(t)+isv(t)ubv(t)=ibv(t)ubv(t)=0。

(8)

5 復(fù)雜動(dòng)態(tài)電路全部支路電流電壓演化規(guī)律求解程序設(shè)計(jì)

依據(jù)前面的討論，我們可以調(diào)用Mathematica強(qiáng)大的符號(hào)矩陣和數(shù)值矩陣運(yùn)算功能，設(shè)計(jì)出Mathematica程序：對(duì)電路中所有信號(hào)源函數(shù)取Laplace變換、列寫(xiě)復(fù)雜多輸入多輸出動(dòng)態(tài)電路s域節(jié)點(diǎn)電壓方程、解出所有支路像電流像電壓函數(shù)、對(duì)支路像電壓像電流函數(shù)取逆Laplace變換給出全部支路電流電壓隨時(shí)間變化的規(guī)律。難解決的任務(wù)都由計(jì)算機(jī)系統(tǒng)來(lái)完成。特別注意，單個(gè)理想電壓源支路串入待定電阻r，適當(dāng)時(shí)間再令r趨于0，是利用節(jié)點(diǎn)電壓法解有理想電壓源支路的復(fù)雜電路的技巧之一。

基于以上分析，我們編寫(xiě)的Mathematica程序主要部分如下：

k=電路支路總數(shù)；n=電路總的節(jié)點(diǎn)數(shù)；T=信號(hào)源周期；

EPTR={{b1,into1,out1,G1,Us1,Is1},{b2,into2,out2,G2,Us2,Is2},……,{bk,intok,outk,Gk,Usk,Isk}}

ivsb=Position[EPTR,1/r][[All,1]];Livsm=Length[ivsb];ivsm={};

For[j=1, j<=Livsm, j++, jb=ivsb[[j]];AppendTo[ivsm,{jb,EPTR[[jb,2]]}]]

Aa=ConstantArray[0,{n,k}];

For[j=1,j

A=Take[Aa,n-1];Clear[Aa] (n-1行k列電路拓?fù)潢P(guān)系矩陣);

Usv=EPTR[[All,5]] (Usv是支路電壓源矢量);

Isv=EPTR[[All,6]] (Isv是支路電流源矢量);

Gdm=DiagonalMatrix[EPTR[[All,4]]] (Gdm是對(duì)角電導(dǎo)矩陣);

AT=Transpose[A] (AT是電路拓?fù)潢P(guān)系矩陣A的轉(zhuǎn)置);

Insv=A.(Gdm.Usv-Isv) (Insv是流入每個(gè)節(jié)點(diǎn)的電流源矢量);

Gm=A.Gdm.AT (Gm是電路節(jié)點(diǎn)電壓方程電導(dǎo)矩陣);

IG=Inverse[Gm] (IG是節(jié)點(diǎn)電壓方程電導(dǎo)矩陣Gm的逆矩陣);

Unv=Limit[IG.Insv,r->0] (Unv是節(jié)點(diǎn)電壓矢量);

Ubv=AT.Unv (Ubv是支路電壓矢量);

For[j=1,j<=Livsm, j++,qb=ivsm[[j,1]];Ubv[[qb]]=Usv[[qb]]];

Uzv=Ubv-Usv (Uzv是電導(dǎo)元件的電壓矢量);

Ibv=Gdm.Uzv+Isv (Ibv是支路電流像函數(shù)矢量);

For[j=1, j<=Livsm, j++,ib=ivsm[[j,1]];ic=ivsm[[j,2]];cin=Drop[Range[1,k],ib];

Ibv[[ib]]=-A[[ic,ib]] (A[[ic,cin]].Ibv[[cin]])] (計(jì)算理想電壓源支路電流);

Izv=Ibv-Isv (Izv是電導(dǎo)元件電流矢量);

unv[t_]=InverseLaplaceTransform[Unv,s,t] (時(shí)域節(jié)點(diǎn)電壓矢量函數(shù));

ubv[t_]=AT.unv[t] (ubv是時(shí)域支路電壓函數(shù));

ibv[t_]=InverseLaplaceTransform[Ibv,s,t] (ibv是時(shí)域支路電流矢量函數(shù));

t1=T;t2=4T;

Table[Plot[ubv[t][[q]],{t,t1,t2},PlotRange->All],{q,1,k}] (支路電壓波形);

Table[Plot[ibv[t][[q]],{t,t1,t2},PlotRange->All],{q,1,k}] (支路電流波形);

usv[t_]=InverseLaplaceTransform[Usv,s,t];isv[t_]=InverseLaplaceTransform[Isv,s,t];

uzv[t_]=ubv[t]-usv[t];izv[t_]=ibv[t]-isv[t] (uzv，izv是電導(dǎo)元件的時(shí)域電壓、電流矢量函數(shù));

p∑[t_]=ibv[t].ubv[t] (時(shí)域各支路總功率函數(shù)，其值應(yīng)該為0)。

6 復(fù)雜動(dòng)態(tài)電路求解實(shí)例

通過(guò)求解16個(gè)元件、13條支路、6個(gè)節(jié)點(diǎn)、2個(gè)電壓源、2個(gè)電流源構(gòu)成的多輸入多輸出復(fù)雜動(dòng)態(tài)電路(見(jiàn)圖2)，展示我們提出的電路解法思想和程序的用法。約定圖中支路參考電流的流入端是各支路電壓的“+”端。電路的索引已經(jīng)標(biāo)于圖中。

圖2 16個(gè)元件、13條擴(kuò)展支路、6個(gè)節(jié)點(diǎn)的多輸入多輸出復(fù)雜動(dòng)態(tài)電路，b=13，n=6

描述4個(gè)信號(hào)源的程序語(yǔ)句為：

ε[t_]=UnitStep[t];is1[t_]=0.01ε[t];is4[t_]=0.02Cos[100t]ε[t];

us7[t_]=15Cos[200t]ε[t];us12[t_]=20Cos[400t]ε[t];T=π/200。

特別在理想電壓源支路中加入虛擬電阻R7=r，于是各支路元件參數(shù)分別為：

{R1,L2,R3,c4,R5,L6,R7,R8,c9,c10,R11,L12,R13}=

{500,0.02,5 000,10-4,2 000,0.01,r,100,10-5,10-6,200,0.02,4 000}。

由圖2見(jiàn)，擴(kuò)展支路數(shù)k=13,節(jié)點(diǎn)數(shù)n=6。特別注意，支路1上的電流源流向與支路1的電流參考流向相反，支路12上參考電流是從電壓源“-”端流入的，因此，電路的電流源矢量isv(t)和電壓源矢量usv(t)為：

isv[t_]={-is1[t],0,0,is4[t],0,0,0,0,0,0,0,0,0};

usv[t_]={0,0,0,0,0,0,us7[t],0,0,0,0,-us12[t],0}。

參考圖2，得到電路元件參數(shù)及拓?fù)潢P(guān)系表示矩陣EPTR：

EPTR={ {1,1,6,1/R1,0,Isv[[1]]},{2,2,6,1/(s×L2),0,0},{3,2,6,1/R3,0,0},

{4,3,6,s×c4,0,Isv[[4]]},{5,4,6,1/R5,0,0},{6,4,6,1/(s×L6),0,0},

{7,5,6,1/R7,Usv[[7]],0},{8,4,5,1/R8,0,0},{9,3,4,s×c9,0,0},

{10,2,3,s×c10,0,0},{11,1,2,1/R11,0,0},

{12,1,3,1/(s×L12),Usv[[12]],0},{13,3,5,1/R13,0,0}}。

運(yùn)行程序后得到的拓?fù)潢P(guān)系矩陣為：

A={ {1,0,0,0,0,0,0,0,0,0,1,1,0},{0,1,1,0,0,0,0,0,0,1,-1,0,0},

{0,0,0,1,0,0,0,0,1,-1,0,-1,1},{0,0,0,0,1,1,0,1,-1,0,0,0,0},

{0,0,0,0,0,0,1,-1,0,0,0,0,-1}}；

ivsm={{7,5}}。

運(yùn)行程序給出電路中支路1—6在時(shí)域[π/200,4π/200]上的電壓電流函數(shù)波形(圖3和圖4)。

圖3 在圖2電路支路1—6的電壓函數(shù)時(shí)域波形

圖4 在圖2電路中支路1—6的電流函數(shù)時(shí)域波形

原電路總的消耗和存儲(chǔ)功率pz(t)的波形和電源提供的總功率ps(t)的波形見(jiàn)圖5，兩個(gè)波形完全相同，電路能量時(shí)時(shí)守恒。

(a)總的消耗和存儲(chǔ)的功率波形 (b)信號(hào)源提供的總功率波形圖5 原電路總的消耗和存儲(chǔ)的功率波形和電源提供的總功率波形

7 結(jié) 語(yǔ)

通過(guò)用Mathematica程序分析方法可以快速解決復(fù)雜系統(tǒng)的動(dòng)態(tài)電路求解問(wèn)題，并且可以直觀形象地以圖像的形式展現(xiàn)其中電流電壓隨時(shí)間變化的關(guān)系。這為生產(chǎn)及科學(xué)研究工作中相關(guān)問(wèn)題的探討提供了一種有效的解決途徑。

[參考文獻(xiàn)]

[1] 王斌,王迎飛,龍姝明.動(dòng)態(tài)電路計(jì)算機(jī)解法研究[J].中國(guó)新技術(shù)新產(chǎn)品,2009(1):10-11.

[2] 張志良.用等效阻抗法解暫態(tài)過(guò)程[J].浙江師大學(xué)報(bào):自然科學(xué)版,2001,24(3):314-316.

[3] 付奎,楊俊海,王風(fēng)華,等.大規(guī)模實(shí)時(shí)電路程序分析方法[J].陜西理工學(xué)院學(xué)報(bào):自然科學(xué)版,2010,26(3):13-17.

[4] 龍姝明,朱杰武,孫彥清,等.數(shù)學(xué)物理方法& Mathematica[M].西安:陜西人民教育出版社,2002:343-367.