《高等数值分析》教学案例的建设与思考

唐玲艳，　宋松和

(国防科学技术大学理学院，湖南长沙410073)



《高等数值分析》教学案例的建设与思考

唐玲艳，宋松和

(国防科学技术大学理学院，湖南长沙410073)

[摘要]《高等数值分析》是一门与实际联系紧密的数学公共课程，它理论深刻，应用广泛.本文结合实际应用，为《高等数值分析》中常微分方程数值解部分设计了一个教学案例，通过理论分析和数值实验向学生展示了刚性问题的概念和相关数值方法，并对《高等数值分析》课程教学案例的设计进行了思考.

[关键词]高等数值分析； 教学案例； 常微分方程数值解； 刚性稳定

1引言

随着计算机技术的飞速发展，科学计算的地位日趋重要，已与实验和理论研究并列，成为现代科学研究的三大主要手段之一.数值分析作为科学计算的基础，渗透到了极其广泛的专业领域.一些重要而基础的数值计算方法，如解线性方程组的Gauss消元法和简单迭代法，解常微分方程组的Euler法、梯形法和Runge-Kutta法等，已为越来越多的工程技术人员所掌握并熟练使用[1-3].与数学家们相比，工程技术人员对这些方法的使用往往更为大胆，这一点在研究生公共课《高等数值分析》的教学过程中也有所反映.理工科研究生普遍对于算法的构造原理，如用有限近似无限，将非线性问题线性化等，接受得很快，编程计算能力也较强，但却容易忽略算法的理论分析，对精度、稳定性和收敛性等概念理解不够深刻.因此在工程应用中遇到一些看似可以算，实则不能算或计算出错的情况时，他们往往无法查找原因并提出解决策略.针对上述问题，我们结合实际应用，精心设计了一个《高等数值分析》的教学案例.一方面向学生揭示蕴含在算法背后的深刻理论，便于其更好地理解和掌握相关知识；另一方面让学生通过研究讨论案例，身临其境地体会实际科研中可能遇到的情况，提高其分析问题和解决问题的能力.

2问题描述：B-Z反应与刚性

一般的化学反应，反应物和产物的浓度单调地发生变化，最终达到不随时间变化的平衡状态.然而在某些反应体系中，有些组分的浓度会忽高忽低，呈现周期性变化，这种现象称为化学震荡.上世纪六十年代初，俄国化学家Belousov和Zhabotinskii首次报道了以金属铈作催化剂，柠檬酸在酸性条件下被溴酸钾氧化时所呈现的化学振荡现象：溶液在无色和淡黄色两种状态间进行着规则的周期振荡，人们将其称为Belousov-Zhabotinskii反应(简称B-Z反应).1972年，Field，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)

为了保证数值解的稳定性，要求‖Yn+1‖≤‖Yn‖.当λ(t)≡λ时，对向前Euler法，即要求-2≤λh≤0.对向后Euler法和梯形法，即要求λh≤0.

可以看出，当λ的绝对值很大时，向前Euler法的步长要取得非常小才能给出满意的结果，向后Euler法和梯形法则没有这样的限制.通常，将像向前Euler法这样稳定区间具有有限边界的方法称为条件稳定，将像向后Euler法和梯形法这样稳定区间为整个负半平面的方法称为A稳定.显然，求解刚性微分方程，最好选用对步长不加限制的A稳定方法.

进一步，考察向后Euler法和梯形法的区别，是否阶数高的方法计算效果就一定好呢？事实并非如此.对变系数情况，向后Euler法和梯形法的稳定性条件分别为

(12)

和

(13)

其中不等式(12)仍未对步长加以限制，不等式(13)在某些条件下相当于

(14)

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

(15)

当λ(t)绝对值很大而又迅速增长时，(15)式是对步长限制非常严苛的条件.这就是说，对于非线性方程情况，梯形法可能要取比向后Euler法小得多的步长，才能获得稳定的结果.

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

假设F(t，Y)关于Y满足Lipschitz条件，Lipschitz常数为L，则

‖G(X)-G(Y)‖=‖F(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，隐式方法的内迭代统一采用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.此时系统的刚性比更大，本文提供的三种方法均会失效，需要寻求更为精细的数值方法.上述分析和计算表明，求解刚性BZ反应方程最好选用对步长不加限制的A稳定，甚至L稳定方法.为了编程和计算的简便，方法最好还是显式的.然而，Dahlquist已经证明：A稳定的方法一定是隐式的，且阶数不超过2[2].为此，Gear放宽A稳定性的要求，提出刚性稳定的概念[5].他证明，当k≤6时，具有σ(μ)=μk形式的k步k阶方法是刚性稳定的，提出著名的Gear方法.前文中的向后Euler法正是一阶Gear方法.此外，由Rosenbrock提出的半隐式Runge-Kutta方法[6]在刚性方程的数值计算中也有很重要的应用.

图2　数值实验的结果，其中ε=0.01，p=6，q=8.4×10-6，v=0.6

6关于教学案例建设的思考

《高等数值分析》是一门与实际联系紧密的数学公共课，它理论深刻，应用广泛.在该课程中采用案例教学，可以最大程度地调动学生学习的积极性，提高其分析问题和解决问题的能力.本文结合实际应用，为《高等数值分析》中常微分方程数值解部分设计了一个教学案例，生动而形象地向学生展示了刚性问题的概念和相关数值方法.结合《高等数值分析》课程的特点，在教学案例的选择和设计过程中，本文主要考虑了以下几点：

( i) 教学案例的选择要大众化，既要有一定的背景，以激发学生的兴趣，又要难度适中，便于学生上手.例如，本文的教学案例选自化学反应中有名的BZ振荡，形式上是一个一阶常微分方程组，课程中介绍的基本方法都可以应用于其上.

( ii) 教学案例一般要结合一定理论，要通过各种信息、知识、经验、观点的碰撞来达到启示理论和启迪思维的目的.例如，本文的教学案例不仅是基本方法在常微分方程组上的运用，也涉及到刚性问题和刚性稳定的概念，学生通过理论分析和数值实验可加深对这部分内容的理解.

( iii) 教学案例的设计要具有层次，既要提供在课程上能够阐述清楚的简单解决方案，又要有一定的可扩展性，能够跟较前沿的研究接轨，指导学生开展课外阅读和进一步的研究.

[参考文献]

[1]徐绪海，朱方生.刚性微分方程的数值解法[M].武汉：武汉大学出版社，1997.

[2]余德浩，汤华中.微分方程数值解法[M].北京：科学出版社，2003.

[3]宋松和，朱建民，成礼智，唐玲艳.高等数值分析课程教学改革探讨[J].高等教育研究学报，2008，31(4) : 66-67.

[4]李如生.非平衡态热力学和耗散结构[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