《高等數(shù)值分析》教學(xué)案例的建設(shè)與思考

唐玲艷，　宋松和

(國(guó)防科學(xué)技術(shù)大學(xué)理學(xué)院，湖南長(zhǎng)沙410073)

?

《高等數(shù)值分析》教學(xué)案例的建設(shè)與思考

唐玲艷，宋松和

(國(guó)防科學(xué)技術(shù)大學(xué)理學(xué)院，湖南長(zhǎng)沙410073)

[摘要]《高等數(shù)值分析》是一門(mén)與實(shí)際聯(lián)系緊密的數(shù)學(xué)公共課程，它理論深刻，應(yīng)用廣泛.本文結(jié)合實(shí)際應(yīng)用，為《高等數(shù)值分析》中常微分方程數(shù)值解部分設(shè)計(jì)了一個(gè)教學(xué)案例，通過(guò)理論分析和數(shù)值實(shí)驗(yàn)向?qū)W生展示了剛性問(wèn)題的概念和相關(guān)數(shù)值方法，并對(duì)《高等數(shù)值分析》課程教學(xué)案例的設(shè)計(jì)進(jìn)行了思考.

[關(guān)鍵詞]高等數(shù)值分析； 教學(xué)案例； 常微分方程數(shù)值解； 剛性穩(wěn)定

1引言

隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展，科學(xué)計(jì)算的地位日趨重要，已與實(shí)驗(yàn)和理論研究并列，成為現(xiàn)代科學(xué)研究的三大主要手段之一.數(shù)值分析作為科學(xué)計(jì)算的基礎(chǔ)，滲透到了極其廣泛的專業(yè)領(lǐng)域.一些重要而基礎(chǔ)的數(shù)值計(jì)算方法，如解線性方程組的Gauss消元法和簡(jiǎn)單迭代法，解常微分方程組的Euler法、梯形法和Runge-Kutta法等，已為越來(lái)越多的工程技術(shù)人員所掌握并熟練使用[1-3].與數(shù)學(xué)家們相比，工程技術(shù)人員對(duì)這些方法的使用往往更為大膽，這一點(diǎn)在研究生公共課《高等數(shù)值分析》的教學(xué)過(guò)程中也有所反映.理工科研究生普遍對(duì)于算法的構(gòu)造原理，如用有限近似無(wú)限，將非線性問(wèn)題線性化等，接受得很快，編程計(jì)算能力也較強(qiáng)，但卻容易忽略算法的理論分析，對(duì)精度、穩(wěn)定性和收斂性等概念理解不夠深刻.因此在工程應(yīng)用中遇到一些看似可以算，實(shí)則不能算或計(jì)算出錯(cuò)的情況時(shí)，他們往往無(wú)法查找原因并提出解決策略.針對(duì)上述問(wèn)題，我們結(jié)合實(shí)際應(yīng)用，精心設(shè)計(jì)了一個(gè)《高等數(shù)值分析》的教學(xué)案例.一方面向?qū)W生揭示蘊(yùn)含在算法背后的深刻理論，便于其更好地理解和掌握相關(guān)知識(shí)；另一方面讓學(xué)生通過(guò)研究討論案例，身臨其境地體會(huì)實(shí)際科研中可能遇到的情況，提高其分析問(wèn)題和解決問(wèn)題的能力.

2問(wèn)題描述：B-Z反應(yīng)與剛性

一般的化學(xué)反應(yīng)，反應(yīng)物和產(chǎn)物的濃度單調(diào)地發(fā)生變化，最終達(dá)到不隨時(shí)間變化的平衡狀態(tài).然而在某些反應(yīng)體系中，有些組分的濃度會(huì)忽高忽低，呈現(xiàn)周期性變化，這種現(xiàn)象稱為化學(xué)震蕩.上世紀(jì)六十年代初，俄國(guó)化學(xué)家Belousov和Zhabotinskii首次報(bào)道了以金屬鈰作催化劑，檸檬酸在酸性條件下被溴酸鉀氧化時(shí)所呈現(xiàn)的化學(xué)振蕩現(xiàn)象：溶液在無(wú)色和淡黃色兩種狀態(tài)間進(jìn)行著規(guī)則的周期振蕩，人們將其稱為Belousov-Zhabotinskii反應(yīng)(簡(jiǎn)稱B-Z反應(yīng)).1972年，F(xiàn)ield，Koros和Noyes用約20個(gè)化學(xué)方程式解釋B-Z反應(yīng)的動(dòng)力學(xué)機(jī)制，提出Oregonator數(shù)學(xué)模型[4]：

(1)

得到無(wú)量綱化后的反應(yīng)動(dòng)力學(xué)方程組

(2)

取y1(0)=4，y2(0)=1.0，y3(0)=4，關(guān)于v的值，Noyes等人建議可取作約0.5.該方程組由三個(gè)非線性微分方程組成，形式復(fù)雜，難以求出解析解，通常借助數(shù)值分析方法對(duì)其進(jìn)行模擬.

3模型分析

為敘述方便，令

則方程組(2)可簡(jiǎn)記為

(3)

Jacobi矩陣為

(4)

(5)

的通解形式為

假設(shè)|Re(λt)|<0 (i=1，2，3)，則t→∞時(shí)，kieλitci→0，Z(t)逐漸收斂到φ(t).并且，|Re(λt)|較大的項(xiàng)衰減快，|Re(λt)|較小的項(xiàng)衰減慢.因此，將

λ1≈-15.611，λ2≈5.0857，λ3≈0.0924.

λ1≈-531.6654，λ2≈26.3259，λ3≈0.0002，

4數(shù)值離散

取分點(diǎn)tn=nh，(n=0，1，2，…)，h為步長(zhǎng)，求解一階常微分方程的常用數(shù)值方法都可用于方程(3)的求解.例如，向前Euler法的計(jì)算公式為

Yn+1=Yn+F(tn，Yn)h.

(6)

向后Euler法的計(jì)算公式為

Yn+1=Yn+F(tn+1，Yn+1)h.

(7)

梯形法的計(jì)算公式為

(8)

以上三種均為單步方法.其中，向前Euler法為一階顯格式，向后Euler法為一階隱格式，梯形法為二階隱格式.假設(shè)F(t，Y)=λ(t)Y，則三者的遞推公式分別為

向前Euler法：

(9)

向后Euler法：

(10)

梯形法：

(11)

為了保證數(shù)值解的穩(wěn)定性，要求‖Yn+1‖≤‖Yn‖.當(dāng)λ(t)≡λ時(shí)，對(duì)向前Euler法，即要求-2≤λh≤0.對(duì)向后Euler法和梯形法，即要求λh≤0.

可以看出，當(dāng)λ的絕對(duì)值很大時(shí)，向前Euler法的步長(zhǎng)要取得非常小才能給出滿意的結(jié)果，向后Euler法和梯形法則沒(méi)有這樣的限制.通常，將像向前Euler法這樣穩(wěn)定區(qū)間具有有限邊界的方法稱為條件穩(wěn)定，將像向后Euler法和梯形法這樣穩(wěn)定區(qū)間為整個(gè)負(fù)半平面的方法稱為A穩(wěn)定.顯然，求解剛性微分方程，最好選用對(duì)步長(zhǎng)不加限制的A穩(wěn)定方法.

進(jìn)一步，考察向后Euler法和梯形法的區(qū)別，是否階數(shù)高的方法計(jì)算效果就一定好呢？事實(shí)并非如此.對(duì)變系數(shù)情況，向后Euler法和梯形法的穩(wěn)定性條件分別為

(12)

和

(13)

其中不等式(12)仍未對(duì)步長(zhǎng)加以限制，不等式(13)在某些條件下相當(dāng)于

(14)

若λ(tn+1)≤λ(tn)，則(14)是成立的；若λ(tn+1)>λ(tn)，則

(15)

當(dāng)λ(t)絕對(duì)值很大而又迅速增長(zhǎng)時(shí)，(15)式是對(duì)步長(zhǎng)限制非常嚴(yán)苛的條件.這就是說(shuō)，對(duì)于非線性方程情況，梯形法可能要取比向后Euler法小得多的步長(zhǎng)，才能獲得穩(wěn)定的結(jié)果.

‖G(X)-G(Y)‖<α‖X-Y‖.

假設(shè)F(t，Y)關(guān)于Y滿足Lipschitz條件，Lipschitz常數(shù)為L(zhǎng)，則

‖G(X)-G(Y)‖=‖F(xiàn)(tn+1，X)-F(tn+1，Y)‖h

(16)

由此推出不動(dòng)點(diǎn)迭代收斂的條件是h<1/L，當(dāng)F(t，Y)關(guān)于Y的變化較劇烈，即L較大時(shí)，條件(16)也是比較嚴(yán)苛的.

5數(shù)值實(shí)驗(yàn)與分析

對(duì)BZ反應(yīng)方程(2)，選取兩組不同的參數(shù)在個(gè)人計(jì)算機(jī)上進(jìn)行數(shù)值實(shí)驗(yàn).機(jī)器型號(hào)為Intel(R)Core(TM)2DuoCPUT5870 2.00GHz2.00GHz，RAM2.86GB，隱式方法的內(nèi)迭代統(tǒng)一采用Newton迭代法，計(jì)算結(jié)果如圖1和圖2所示.可以看出：

圖1　數(shù)值實(shí)驗(yàn)的結(jié)果，其中ε=0.03，p=2，q=0.006，h=0.5

(iii) 事實(shí)上，化學(xué)測(cè)量的結(jié)果顯示[1]，BZ反應(yīng)方程(2)中各參數(shù)的量級(jí)約為ε～2×10-4，p～3.1×102，q～8.4×10-6，v～0.5.此時(shí)系統(tǒng)的剛性比更大，本文提供的三種方法均會(huì)失效，需要尋求更為精細(xì)的數(shù)值方法.上述分析和計(jì)算表明，求解剛性BZ反應(yīng)方程最好選用對(duì)步長(zhǎng)不加限制的A穩(wěn)定，甚至L穩(wěn)定方法.為了編程和計(jì)算的簡(jiǎn)便，方法最好還是顯式的.然而，Dahlquist已經(jīng)證明：A穩(wěn)定的方法一定是隱式的，且階數(shù)不超過(guò)2[2].為此，Gear放寬A穩(wěn)定性的要求，提出剛性穩(wěn)定的概念[5].他證明，當(dāng)k≤6時(shí)，具有σ(μ)=μk形式的k步k階方法是剛性穩(wěn)定的，提出著名的Gear方法.前文中的向后Euler法正是一階Gear方法.此外，由Rosenbrock提出的半隱式Runge-Kutta方法[6]在剛性方程的數(shù)值計(jì)算中也有很重要的應(yīng)用.

圖2　數(shù)值實(shí)驗(yàn)的結(jié)果，其中ε=0.01，p=6，q=8.4×10-6，v=0.6

6關(guān)于教學(xué)案例建設(shè)的思考

《高等數(shù)值分析》是一門(mén)與實(shí)際聯(lián)系緊密的數(shù)學(xué)公共課，它理論深刻，應(yīng)用廣泛.在該課程中采用案例教學(xué)，可以最大程度地調(diào)動(dòng)學(xué)生學(xué)習(xí)的積極性，提高其分析問(wèn)題和解決問(wèn)題的能力.本文結(jié)合實(shí)際應(yīng)用，為《高等數(shù)值分析》中常微分方程數(shù)值解部分設(shè)計(jì)了一個(gè)教學(xué)案例，生動(dòng)而形象地向?qū)W生展示了剛性問(wèn)題的概念和相關(guān)數(shù)值方法.結(jié)合《高等數(shù)值分析》課程的特點(diǎn)，在教學(xué)案例的選擇和設(shè)計(jì)過(guò)程中，本文主要考慮了以下幾點(diǎn)：

( i) 教學(xué)案例的選擇要大眾化，既要有一定的背景，以激發(fā)學(xué)生的興趣，又要難度適中，便于學(xué)生上手.例如，本文的教學(xué)案例選自化學(xué)反應(yīng)中有名的BZ振蕩，形式上是一個(gè)一階常微分方程組，課程中介紹的基本方法都可以應(yīng)用于其上.

( ii) 教學(xué)案例一般要結(jié)合一定理論，要通過(guò)各種信息、知識(shí)、經(jīng)驗(yàn)、觀點(diǎn)的碰撞來(lái)達(dá)到啟示理論和啟迪思維的目的.例如，本文的教學(xué)案例不僅是基本方法在常微分方程組上的運(yùn)用，也涉及到剛性問(wèn)題和剛性穩(wěn)定的概念，學(xué)生通過(guò)理論分析和數(shù)值實(shí)驗(yàn)可加深對(duì)這部分內(nèi)容的理解.

( iii) 教學(xué)案例的設(shè)計(jì)要具有層次，既要提供在課程上能夠闡述清楚的簡(jiǎn)單解決方案，又要有一定的可擴(kuò)展性，能夠跟較前沿的研究接軌，指導(dǎo)學(xué)生開(kāi)展課外閱讀和進(jìn)一步的研究.

[參考文獻(xiàn)]

[1]徐緒海，朱方生.剛性微分方程的數(shù)值解法[M].武漢：武漢大學(xué)出版社，1997.

[2]余德浩，湯華中.微分方程數(shù)值解法[M].北京：科學(xué)出版社，2003.

[3]宋松和，朱建民，成禮智，唐玲艷.高等數(shù)值分析課程教學(xué)改革探討[J].高等教育研究學(xué)報(bào)，2008，31(4) : 66-67.

[4]李如生.非平衡態(tài)熱力學(xué)和耗散結(jié)構(gòu)[M].北京：清華大學(xué)出版社，1986.

[5]GearC.William.Theautomaticintegrationofstiffordinarydifferentialequations[J].InformationProcessing68，ed.，Morrell，AJH:NorthHallandPublishingCo.， 1969: 187-193.

[6]AnFeng，CharlesD.Holland，StevenE.Gallun.Developmentandcomparisonofageneralizedsemi-implicitRunge-KuttamethodwithGear’smethodforsystemsofcoupleddifferentialandalgebraicequations[J].Computer&ChemicalEngineering， 1984， 8(1): 51-59.

ConstructionandThinkingofTeachingCasesof

AdvancedNumericalAnalysis

TANG Ling-yan，SONG Song-he

(DepartmentofMathematicsandSystemScience，ScienceSchool，NationalUniversityof

DefenceTechnology，Changsha410073，China)

Abstract:Advancednumericalanalysisisapublicmathematicalcoursewhichislinkedcloselywithpractice.Ithasprofoundtheoryandwideapplications.Combinedwiththepracticalapplication，ateachingcaseofthenumericalsolutionofordinarydifferentialequationsisdesignedforadvancednumericalanalysisinthispaper.Throughtheoreticalanalysisandnumericalexperimentsoftheteachingcase，weshowstudentstheconceptofstiffproblemsandrelatednumericalmethods.Attheendofthepaper，somethinkingismadefortheconstructionofteachingcasesofadvancednumericalanalysis.

Keywords:advancednumericalanalysis;teachingcases;numericalsolutionofordinarydifferentialequations;stiffstability

[基金項(xiàng)目]研究生數(shù)學(xué)公共課一流課程體系建設(shè)項(xiàng)目; 國(guó)防科學(xué)技術(shù)大學(xué)校本科教改項(xiàng)目

[收稿日期]2014-11-12

[中圖分類號(hào)]G642.0

[文獻(xiàn)標(biāo)識(shí)碼]C

[文章編號(hào)]1672-1454(2015)01-0042-06