求解指標(biāo)1的微分代數(shù)方程組的一類新方法

鮑文娣, 韓海力

(1.中國石油大學(xué) 理學(xué)院, 山東 青島 266552； 2.青島杰瑞自動(dòng)化有限化司, 山東 青島 266071)

求解指標(biāo)1的微分代數(shù)方程組的一類新方法

鮑文娣1, 韓海力2

(1.中國石油大學(xué) 理學(xué)院, 山東 青島 266552； 2.青島杰瑞自動(dòng)化有限化司, 山東 青島 266071)

給出求解指標(biāo)1的微分代數(shù)方程組的一類新的計(jì)算方法．將微磁學(xué)仿真的方法推廣到求解微分代數(shù)方程組,并給出方法的收斂性和相容性分析． 利用與伴隨法相復(fù)合的方法,提高方法的收斂階．并將方法應(yīng)用于晶體管放大器的模型中．?dāng)?shù)值實(shí)驗(yàn)表明方法是有效的．

Gauss-Seidel方法; 微分代數(shù)方程組; 伴隨法; 投影法

0 引言

許多重要的數(shù)學(xué)模型都表示為微分代數(shù)方程組的形式, 它在各式各樣的科學(xué)工程都有應(yīng)用．例如:大規(guī)模電路分析、計(jì)算機(jī)輔助設(shè)計(jì)、多體動(dòng)力系統(tǒng)的實(shí)時(shí)仿真、能源系統(tǒng)、化學(xué)反應(yīng)仿真、最優(yōu)控制等等．很多物理中的問題最初的模型即為一個(gè)微分代數(shù)方程組的形式．

近些年來,微分代數(shù)方程組的求解已成為近十年來研究的熱點(diǎn)問題．一些求解微分代數(shù)方程組的數(shù)值方法得到了進(jìn)一步的發(fā)展,如Padé逼近、波行迭代法[2-10]等.

本文提出一類新的數(shù)值方法,探討應(yīng)用復(fù)合的Gauss-Seidel法和SOR方法來求解微分代數(shù)方程組,并進(jìn)一步給出方法的收斂性分析．在第2部分構(gòu)造適用于微分代數(shù)方程組的復(fù)合Gauss-Seidel法和SOR方法,并給出了方法的誤差估計(jì)．在第三節(jié)中,就用于求解晶體管放大器模型和擺模型中．最后,在第四節(jié)中給出結(jié)論.

1 求解微分代數(shù)方程的數(shù)值方法

1.1 非線性微分代數(shù)方程組

考慮下面的微分代數(shù)方程組:

(1)

其中x∈Rnx,y∈Rny,f∈Rnx,g∈Rny,并且設(shè)(1)中的微分代數(shù)方程組指標(biāo)為1,也就是說g(x(t),y(t))關(guān)于y的Jacobian矩陣gy在(1)的解的鄰域內(nèi)是可逆的．

眾所周知,一般的前Euler方法(2)求解微分代數(shù)方程組時(shí),是不穩(wěn)定的．下面我們用考慮Gauss-Seidel方法來修正(2)式．

1) 前Euler方法:

(2)

2) 顯式的Gauss-Seidel方法:

(3)

3) 隱式的Gauss-Seidel方法:

(4)

4) SOR方法

(5)

通過比較發(fā)現(xiàn)上述四種方法中,隱式的方法最穩(wěn)定．為了提高誤差階,將復(fù)合法來求解指標(biāo)1的微分代數(shù)方程組．利用Taylor展開式,可得

yn-y(tn)=O(h),xn-x(tn)=O(h),

yn-y(tn)=O(h2),xn-x(tn)=O(h2)．

若Jacobian矩陣gy滿足Lipschitz條件‖gy(x,φ(x))‖≤L,則由文[10]中的定理3.4可知全局誤差E=x(t)-xN滿足:

(6)

其中t0為初始時(shí)間,xN是第N步的數(shù)值解,C是獨(dú)立于h的常數(shù)．(6)式說明復(fù)合法是一階精度的全局收斂的方法．

1.2 線性微分代數(shù)方程組

考慮下面的微分代數(shù)方程組

(7)

矩陣B(t)分裂成B(t)=D(t)-L(t)-U(t),其中D(t)是對角元素,L(t)和U(t)分別為下三角和上三角矩陣．利用Euler方法,(7)式可寫成下面的形式:

由SOR迭代,有

(h(D-wL)+wA)x(k+1)=(h(wU+(1-w)D)+wA)xk+whf

(8)

其中w為參數(shù)．若h(D-wL)+wA可逆,則問題有唯一的數(shù)值解x(k+1),并且只有當(dāng)(wA,D(t)-wL(t))是正則時(shí)才有可能．所以為了方法能正常使用,需要假設(shè)矩陣對(wA,D(t)-wL(t))對于所有的t∈I是正則的．

注1 1) 當(dāng)w=1時(shí),(8)為Gauss-Seidel法;2) 方法(7)的誤差估計(jì)和穩(wěn)定性可由Taylor展開得到,經(jīng)計(jì)算局部誤差為O(h2)；3) 對于大規(guī)模方程組方法(8)簡單有效．

2 數(shù)值實(shí)例

例1 晶體管放大器(圖1),圖中Ue(t)為輸入電壓,Ub(t)為工作電壓,Ui(t)(i=1,2,3,4,5)分別為節(jié)點(diǎn)1,2,3,4,5處的電壓,U5(t)為輸出電壓．電流通過電阻時(shí)滿足I=U/R, 電流通過電容時(shí)有I=C·dU/dt,其中R,C為常數(shù),U為電壓．晶體管是這樣工作的,將來自節(jié)點(diǎn)3、4 的電流放大到節(jié)點(diǎn)2到節(jié)點(diǎn)3的電流的99倍,這依賴于電壓差U3-U2．Kirchhoff表述為在任一時(shí)刻匯集于電路中任一節(jié)點(diǎn)的各支路電流的代數(shù)和等于零．將這一法則應(yīng)用于圖1中的節(jié)點(diǎn)5可得下面的方程組:

(9)

初始信號(hào)為:Ue=0.4sin(200πt)．由穩(wěn)式的Gauss-Seidel 法復(fù)合而得到的復(fù)合法應(yīng)用到晶體放大器問題(9)中,其數(shù)值解由圖2給出,比較輸入U(xiǎn)e和輸出電壓U5(t)可知我們的方法是有效的．

圖1 晶體管放大器原理圖

例2 圖2中(擺)懸掛在可忽略質(zhì)量的長為l的棒上的質(zhì)點(diǎn)m,在重力g作用下進(jìn)行擺動(dòng)．設(shè)質(zhì)點(diǎn)所在位置坐標(biāo)為(p,q),則運(yùn)動(dòng)滿足下面的方程:

(10)

(11)

0=p2+q2-l2

(12)

其中(u,v)是速度,λ是棒的張力．表達(dá)式(6)是指標(biāo)為3的微分代數(shù)方程組．微分一次可得:

0=pu+qv

(13)

從幾何上理解為速度為流形(13)的切線,也就是說正交于梯度2(p,q)．方程組(10) (11) (13)是指標(biāo)為2的形式．再微分一次并且由(12)可得:

0=m(u2+v2)-gq-l2λ

(14)

則方程組(10) (11) (14)為指標(biāo)1的形式．將顯式的Gauss-Seidel 方法應(yīng)用于帶有下面相容初始條件的微分代數(shù)方程組(10) (11) (14)中,

p(0)=1,q(0)=0,u(0)=0,v(0)=0,λ(0)=0,

圖2 輸入電壓Ue與輸出電壓U5(t)關(guān)系

則有下列關(guān)系式:

(15)

(16)

(17)

顯然公式(15)-(17)是精度為1階的復(fù)合法(文[11]中的可分的Runge-Kutta法,當(dāng)s=1時(shí))．由圖3和圖4可以看出我們的復(fù)合法能得到與文[12]中RADAU5類似的結(jié)果．

圖3 p和它的導(dǎo)數(shù)p′的數(shù)值解

圖4 q和它的導(dǎo)數(shù)q′的數(shù)值解

3 結(jié)論和展望

文中給出了求解微分代數(shù)方程的Gauss-Seidel 方法和 SOR方法．將復(fù)合的Gauss-Seidel方法應(yīng)用到晶體管放大器的模型中,結(jié)果表明方法是有效的．另一方面,將Gauss-Seidel法與投影法相復(fù)合用于求解擺振動(dòng)問題．兩種方法簡單易實(shí)現(xiàn),所以適用于工程計(jì)算中．然而,方法的精度較低,且只用于求解指標(biāo)1的微分代數(shù)方程組．考慮高精度的適用于高指標(biāo)的微分代數(shù)方程組的計(jì)算方法是我們下一步的工作．

[1] Wang X P, Carlos J. Garcia-Cerververa, Weina E. A Gauss-Seidel Projection method for miscromagnetics simulations[J]. Journal of computational physics, 2001, 171:357-372.

[2] Jiang Y, Wing O. Convergence conditions of waveform relaxation methods for circuit simulation[J]. Proceedings of the 1998 IEEE International Symposium on Circuits and Systems, 1998, 6: 232-235.

[3] Leimkuhler B, Miekkala U, Nevanlinna O. Waveform relaxation for linear RC-circuits[J], IMPACT of Computing in Science and Engineering. 1991, 3: 123-145.

[4] Miekkala U. Dynamic iteration methods applied to linear DAE systems[J]. J Comput. Appl.Math. 1989, 25: 133-151.

[5] Song Y. Waveform relaxation methods and accuracy increase[J]. Ann of Diff Eqs, 1995,11: 440-454.

[6] Roberto Barrio. Performance of the Taylor series method for ODEs/DAEs[J]. Applied Mathematics and Computation, 2005,163: 525-545.

[7] Wang H S. Numerical solutions of Hamiltonian DAEs and applications[J]. Applied Mathematics and Computation, 2006, 180: 97-110.

[8] Brenan K E, Campbell S L, Petzold L R. Numerical Solution of Initial Value Problems in Differential Algebraic Equations[M]. Elsevier, New York, 1989.

[9] Hairer E, Norsett S P, Wanner G.Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems[M]. Springer, New York, 1992.

[10] Hairer E, Norsett S P, Wanner G. Solving Ordinary Differential Equations I. Stiff and Differential-Algebraic Problems[M]. Springer, New York, 1992.

[責(zé)任編輯：李春紅]

ANewClassofMethodsforDifferential-algebraicEquationsofIndex-one

BAO Wen-di1, HAN Hai-li2

(1.School of science, China University of Petroleum, Qingdao Shandong 266552, China)(2.Qingdao JARI Automation Co., Qingdao Shandong 266071, China)

In this paper, a new class of methods for solving the index-1 differential-algebraic equations are proposed. The methods are the extensions of gauss-seidel method to DAEs. In order to increase the convergence order, the composition with the adjoint method is provided and then the convergence of the composition method is given. The numerical tests applied to problems from an amplifier and a pendulum show the algorithms are feasible.

gauss-seidel method; differential-algebraic equations (DAEs); adjoint method; projection method

O175.22

A

1671-6876(2012)02-0117-05

2012-02-16

鮑文娣(1979-), 女, 河北武邑人, 講師, 博士研究生, 研究方向?yàn)榫仃囉?jì)算及微分代數(shù)方程組求解.