唐玲艷，　宋松和

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

?

《高等數值分析》教學案例的建設與思考

唐玲艷，宋松和

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

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

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

1引言

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

2問題描述：B-Z反應與剛性

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

(1)

得到無量綱化后的反應動力學方程組

(2)

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

3模型分析

為敘述方便，令

則方程組(2)可簡記為

(3)

Jacobi矩陣為

(4)

(5)

的通解形式為

假設|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數值離散

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

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

(6)

向后Euler法的計算公式為

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

(7)

梯形法的計算公式為

(8)

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

向前Euler法：

(9)

向后Euler法：

(10)

梯形法：

(11)

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

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

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

(12)

和

(13)

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

(14)

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

(15)

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

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

假設F(t，Y)關于Y滿足Lipschitz條件，Lipschitz常數為L，則

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

(16)

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

5數值實驗與分析

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

圖1　數值實驗的結果，其中ε=0.03，p=2，q=0.006，h=0.5

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

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

6關于教學案例建設的思考

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

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

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

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

[參考文獻]

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

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

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

[4]李如生.非平衡態(tài)熱力學和耗散結構[M].北京：清華大學出版社，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

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

[收稿日期]2014-11-12

[中圖分類號]G642.0

[文獻標識碼]C

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