唐玲艷，　宋松和

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

?

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

唐玲艷，宋松和

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

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

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

1引言

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

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

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

(1)

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

(2)

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

3模型分析

為敘述方便，令

則方程組(2)可簡記為

(3)

Jacobi矩陣為

(4)

(5)

的通解形式為

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

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

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

4數(shù)值離散

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

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

(6)

向后Euler法的計算公式為

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

(7)

梯形法的計算公式為

(8)

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

向前Euler法：

(9)

向后Euler法：

(10)

梯形法：

(11)

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

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

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

(12)

和

(13)

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

(14)

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

(15)

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

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

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

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

(16)

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

5數(shù)值實驗與分析

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

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

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

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

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

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

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

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

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

[參考文獻]

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

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

[3]宋松和，朱建民，成禮智，唐玲艷.高等數(shù)值分析課程教學(xué)改革探討[J].高等教育研究學(xué)報，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

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

[收稿日期]2014-11-12

[中圖分類號]G642.0

[文獻標識碼]C

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