（应用数学专业论文）计算机代数在微分方程近似解及其化简中的应用与研究.pdf

摘要 应用过差塑岱兹系统求数学或工程问题中的符号解，进行公式推导是 计算机代数研究的一个重要方面。本文应用m a t h e m a t i e a 系统的强大的 符号运算功能，作图功能，模式匹配功能，以及m a t h e m a t i c a 系统提供的 控制语句，对一类弱非线性系统 6 + 5 * u e * f ( u ，血) 卜 的有效渐近展开式解进行了研究，用m a t h e m a t i e a 系统实现了自动求解 问题。其中所涉及到的方。法是一些古典的奇异摄动方法的求解，如l i n d s t e d t - - p o i n c a r e 方法、l i e 变换方法、多尺度法、平均化方法等，这里大帮询 题只考虑到自治系统的自动求解。有时为了研究的方便，往往要把高维的 微分系统进行化筒变成低维的微分系统，本文的最后一章论述了用m a t h e m a t i c a 系统实现中心流形方法自动化简微分方程的基本思想及过程。最 后调试通过所有程序并做成程序包。该程序包具有良好的扩充性；可以作 为奇异摄动方法的教堂筮鲑。 a b s t ra c t a p p l y i n gc o m p u t e ra l g e b r a i cs y s t e m t og e tf o r m a ls o l u t i o ni np r o b l e mo fm a t h e m a t i c so re n g i n e e r i n ga n df o r m a ld e d u c t i o n sa r ei m p o r t a n t a s p e c t si ns t u d y i n gc o m p u t e ra l g e h r a t h i sp a p e ru s e sp o w e r f u ls y m b o l i co p e r a t i o n ，d r a w i n g ，p a t t e r nm a t c h i n ga n dc o n t r o l s e n t e n c ep r o v i d e db y m a t h e m a t i c at ow e a k l yn o n l i n e a rs y s t e m i i + 碗* u e * f ( u ，d ) ， s t u d i e si t ss o l u t i o n so fe f f e c t i v e l ya s y m p t o t i ce x p a n s i o n s , a n da u t o m a t i c a l l yg e t si t ss o l u t i o n so fa s y m p t o t i ce x p a n s i o n s , w h e r e s o m em e t h o d so f c l a s s i c a l l ys i n g u l a rp e r t u r b a t i o n sa r eu s e di n t h ep a p e r , s u c ha sl i n d s t e d t p o i n c a r e ，l i et r a n s f o r m a t i o n ，m u l t i p l es c a l e s ，a v e r e g i n g ，a n dw h e r e s o m eo fa u t o m a t i c a l l ys o l v i n gp r o b l e m so fa u t o n o m o u ss y s t e m sa r ec o n s i d e r e d s o m e t i m e si t7 sc o n v e r t i e n tt os t u d yr d i n a r yd i f f e r e n t i a le q u a t i o n s t h a tu s u a l l yt r a n s f o r m i n gs y s t e m so fo r d i n a r yd i f f e r e n t i a le q u a t i o n so f b i g h e rd i m e n s i o n si n t os y s t e m so fl o w e rd i m e n s i o n s t h el a s tc h a p t e r d i s c u s s e st h eb a s i ci d e aa n do p e r a t i o n a lp r o c e s so fa p p l y i n gc e n t r a lm a n i f o l dt oa u t o m a t i c a l l yr e d u c i n gd i f f e r e n t i a le q u a t i o n sw i t hm a t h e m a t i c a a tl a s tc o m p i l e dp r o g r a m sa r ei n t e g r a t e di n t op a c k a g e s w h i c hc a nb e e x p a n d e d ，a n dc a nb ec o n s i d e r e da st e a c h i n gs o f t w a r e a b o u t p e r t u r b a t i o n m e t h o d s 中南工业大学硕士学位论文 第一章基本知识与理论 本章阐述了计算机代数的研究对象，意义及发展概；m a t h e m a t i c a 系统软 件的特点及主要功能与应用 筒述了摄动方法求微分方程渐近解的基本思想。 1 1 引亩 1 1 1 计算机代数爰其应用与发展 1 ) 计算机代数的研究对客 计算机代数( c o m p u t e ra l g e b r a ) 又称计算机公式处理或符号处理，是一门 专门研究使用计算机进行公式推演的学科，研究公式推演的算法和计算机语言 与系统。 与数值计算不同数值计算的目的在于更有效地使用计算机，以求得各类 数学问题的满足一定精度的数值解，而计算机代数是由计算机自动推演数学公 式的新技术，是人工智能的一个应用领域，利用计算机代数系统可以获得对公 式处理的结果或各类问题的解析解。 计算机代数所处理的对泉是数学公式，如t 数学表达式的化简，多项式的因 式分解，求两个多项式的最大公因式，符号地进行矩阵的代数运算，求代数方程 的解公式，求微分与积分，求微分方程与积分方程的解析解，数学定理的证明 等。这一点与数值计算有明显的区别，如：a + l 4 = 7 1 2 ，1 2 6 4 = 3 1 6 ，( a + b ) 2 一a b = a 2 + a b + b 2 ，这些是数值计算方法不能得到的 计算机代数所涉及到的知识也非常广泛，如数学中的许多分支，特别是数 论与近世代数等高散数学领域。蔷掌握程序设计语言、数据结构、编译原理、算 法设计与分析等计算机学科知识。它所处理的数学问题不限于代数，还包括微 积分及其它分支由计算机代数的成果使得计算机能自动地进行符号计算，这 样计算机代数便成了研究数学的工具，同时又是以数学为工具的学科，成为其 他学科的研究工具。计算机代数可以用糖确的算法求解数学问题，如用计算机 求原函数，这也是人工智能专家系统的典型成果之一。计算机代数是计算机科 1 中南工业大学硕士学位论文 学、数学、人工智能三者之间的边缘学科。 2 ) 计算机代数的应用 计算机代数可以应用到许多领域，如：电路理论，自动控制系统，定理的规 器证明，人工智能、生物、化学、天体力学、广义相对论、工程设计、教育学等多个 领域。许多系统解决了大量的应用问题，许多系统又在解决应用课题的过程中 得到完善与提高。从射影几何那样的纯粹数学到机器人的手臂运动那样的技术 科学，从化学元素的微量分析到天体力学中的星球运动的宏观计算等，计算机 代数都傲出了杰出的贡献。 目前，世界上已有6 0 多种计算机代数系统，而且功能越来越完善，可靠性越 来越高，例如：前苏联m i r - - 2 型计算机上的a n a l i t i k 系统发现了格拉斯坦 ( g r a d s h t e y n ) 和吕齐克( r y z h i k ) 积分表( 前苏联较新的版本) 的若干印刷 错误。在物理方面。为适应高能物理研究的需要，人们在系统中增加一些软件 包，以扩大狄拉克代数，模型匹配和多项式的处理能力，其中微扰级数的展开， 费曼积分的计算和狄拉克矩阵的处理等大量的公式推导和计算的各个步骤都 可由计算机来实现，面且用手工处理，除非极简单的情况，几乎是不可能的。在 天文学方面，1 8 6 7 年，狄劳奈( d e l a u n a y ) 发表了关于月球轨道计算公式，这个结 果花费了他2 0 多年的心血。最近狄普里特( d e p r i t ) 等人利用罗姆( r o m ) 的计算 机代数系统在小型机上进行复算，只花了2 0 小时，结果表明狄劳奈的结果在第 七级有一个小错误。用计算机代数的方法可以把科学工作者从繁重的数学计算 中解救出来。计算机代数还可以用于生物，化学中分子结构研究，及广义相对论 中的张量的运算等 3 ) 计算机代数的发展 用计算机进行公式处理最早可以追溯到1 9 5 3 年，凯利曼安( h g k a h r i m a n i a n ) 和诺兰( j n o l a n ) 首先使用数字计算机进行解析运算。6 0 年代初， 麻省理工学院的赛( s i n g ) 和邢( s 珊) 在计算机上做了不定积分的实验。1 9 6 9 年 以后，滚边睾郎在日本以常徽分方程为对象的研究取得了新的成果。 自1 9 7 1 年开始。每五年召开一次国际会议，专门交流计算机代数方面的成 果。其间在欧洲还另行召开了两次国际会议。当时的公式处理都是在高速的，大 贮存容量的大型机上运行，1 9 7 8 年默马思( m u m a t h ) 使用超小型的微型计算机 进行公式处理，取得了巨大的成功。以后许多人使用小型机或微型机进行公式 处理均取得了成果，这样计算机代数迅速地实用化、大众化。究其原因( 1 ) 计算 2 中南工业大学硕士学位论文 机的速度与容量大幅度地提高，( 2 ) 高效率的算法的研究与开发以及经验的结 果。这样，计算机代数作为计算机科学中的新分支有待于进一步开发，特别是国 内还落后子国际水平。 1 1 2 文献综述 随着计算机的速度，容量的大幅度地提高，高效率算法不断出现，计算机代 数得到了蓬勃地发展。到目前为止，世界上已有6 0 多种符号计算机系统，其中较 为出名的有m a c s y m a ，r e d u c e ，s c r a t c h p a d ，s a c ，m a t h e m a t i c a 。在国 内，也做了一些计算机代数系统的研制工作，如：中国科学院陈东岳等人研制的 f c y 系统，江西工业大学衷仁保等研刳的c a s c 系统。 计算机代数系统一般具有如下基本功能：多项式和有理数的加、减、乘、除、 乘方等运算，单变量和多变量多项式在整数环上的因式发解，表达式的化筒；矩 阵的运算，初等函数的符号积分与微分各种形式的代换和模式匹配f 代数方程 的求解等。而且大多数系统都是交互式的，用起来直观方便。 计算机代数系统一般可以分为两种。一种是程序设计武系统，如。麻省理工 学院的m a c s y m a 系统f 兰镱公司与犹地大学的r e d u c e 系统。这类系统具 有较强的通用性，可以用来设计一些较好的程序包l 男一类系统本身就好象一 个大的程序包，其实就是一个专家系统，这类系统是专为一些具体问题而设计 的，如a l t r a n 是专门用于有理方程的求解。 m a t h e m a t i c a 系统是美国物理学家s t e p h e n ，w o l f r a m 领导的一个小组开发 的。这是一个集成化的计算机软件系统，茚可以作为专家系统，又可以作用通用 系统。主要功能有三个方面t 符号演算，数值计算，图形。可以求多项武方程，有 理式方程和超越方程的精确解与近似解，还可以求援限，导函数；求幂级数展开 武 求解某些微分方程等功能，可以做长整数，长整数分母的有理效的精确运 算，任意精度的数值运算，并且有强大的数学作图功能。系统本身有很多程序 包，可以求解许多特殊的数学问题因此，m a t h e m a t i c a 系统便成了数学研究与 工程研究中的非常有力的工具。 微分方程的求解是一个很古老的数学问题，雨且有非常好的堡论与算法求 一个微分方程的数值解。若一个微分方程存在闭形式的解析解，又怎样用计算 机求出其解析解呢? 又怎样用计算机求微分方程的近似分析解呢? 常见的方法 是付立叶( f o u r i e r ) 变换与控普拉斯( l a p l a c e ) 变换( 见文献e s e e 3 s ) ，积分等方 3 中南工业大学硕士学位论文 法。这些方法一般很难在计算机上实现。由于m a t h e m a t i c a 系统有很强的符号 推演能力，级数运算能力，这样就可以方便地对一类特殊的常微方程进行自动 地求近似级数解。 在解决工程技术和科学领域中的各种理论与实际问题时，就数学方面来 说，主要是求分析解与数值解。由于数字计算机的出现，数值计算得到了飞速地 发展，但这并不能取代近似分析解，摄动方法可以得到简单而又准确的近似分 析解，且容易看出每个参数对解的影响，弄清解的结构，由于计算机代数系统有 很强的符号计算能力，这样计算机代数便成了摄动方法中的强有力工具，基于 m a t h e m a t i c a 的符号运算的功能，可以方便地求出一类弱非线系统的摄动近似 解。 摄动方法是近百年来由许多数学家、力学家及物理学家共同建立和发展起 来的具有广泛应用的一个重要方法。最常见的奇异摄动方法有：应变坐标法、匹 配法、平均法、多重尺度法。在十九世纪天文学家l i n d s t e d t 等为了求方程 1 i + ( 1 】l u = f ( u ，矗)e 一( 丢) p l o t c 表达式( 或表达式表) ， 交量，上限，下限) ，可选项 ，按变量所定义的 范围，作出表达式的数学图形。如t i f l 8 ；= p l o t s i n x * c o s x ，( x ，一p i ，p i ) o u t 8 = 小 0 4 o 2 ， 一3 2 一1 1 v - o u 4 中南工业大学硕士学位论文 一g r a p h i c s m d r o p 表，整数n 或d r o p f 表， 整数m ，整数n ) ，从原表中删除元素。 p r e p e n d - 表，表达式 ，把表达式放在原表所有元素的前面，如； i n 9 ：= p r e p e n d f ( 1 ，2 ，3 ) ，o o u t f 9 = 0 ，l ，2 ，3 f l a t t e n - 表 ，把原表抹平得到一层表。 f r e e q 表达式，子表达式 ，判断子表达式是不是在表达式中，不在则返回 真，否则返回假。 m a p 函数，表达式 ，用函数作用在表达式的第一层的每一个元素上，由作 用的结果构成新的表达式，如： i n 1 0 t = m a p e s i n ， l ，2 ，3 ) 3 o u t 1 0 = ( s t u d ，s m e 2 ，s i n 3 ) 2 ) 控制语句与模式，变换规则 m a t h e m a t i c a 系统对表达式的处理过程一般统称为求值。求值的对象是表 达武，最终得到结果也是表达式，求值是一个表达式到另一个表达式的映射。这 样在m a t h e m a t i c a 系统内部就有一套非常丰富的映射规则，即内部规则库。综上 所述，m a t h e m a t i c a 系统的表达武处理系统是由一个求值系统和一个变换规则 库组成，见下图： 值系统 m a t h e m a t i e a 的求值过程 为了编程的方便，m a t h e m a t i c a 系统还允许用户自己定义规则库，这便于定 义用户自己的函数，如定义函数f 8 中南工业大学硕士学位论文 f x 一 ：= x + 3 那么当x 赋值为5 时，便有 f e s 3 = s 其中“x 一”中的下画线在系统里称为空白，它是m a t h e m a t i c a 系统中最基本的模 式，模式在m a t h e m a t i c a 系统中的作用非常重要，它提供了m a t h c m a t i c a 系统中 有关表达式匹配的原则，如s i n 能与任意正弦函数匹配，即s i n x 、s i n l - 3 3 等 匹配。 m a t h e m a t i c a 系统的规则分为两类，自动使用规则与非自动使用规则。“= ” 与“。= ”都孀于描述变换规则定义的符号，都是用于定义自动使用规则。非自动 规则的描述见如下； 模式- - 表达式右部表达式在定义时求值 模式； 表达武右部表达式在定义时不求值。 m a t h c m a t i c a 系统也象其他程序设计语言一样提供控制语句。m a t h e m a t i c a 程序的执行分为三种；顺序执行、选择执行、循环执行，即提供了i f 语句，f o r ， w h i l e ，d o 循环语句。也提供c o n t i n u e ，b r e a k ，r e t u r n 等语句。同时m a t h e m a t i - c a 也提供引进局部变量的机制，可以使用m o d u l e ( b l o c k ) 结构定义局部变量， m o d u l e 结构的使用形式是： m o d u l e e 局部变量表，表达式 。 当不要局部变量时，可以使用纯函数形式。 1 2 3m a t h e m a t i c a 系统在工程数学中的举例 作为例子。考虑一个阶常系数常微分方程组： f x i 一缸z + ( 1 ，2 ，1 ) i x 2 = 3 x 1 + 3 x 2 在m a t h e m a t i c a 系统里面有几种方法可以求解该方程组，首先把常微分系统写 成一表的形式。并赋给一变量t e m p ： l n l l ：= t e m p = x 1 7 t ；- - - - 4 * x l t + 2 * x l t 。 x 2 , e t 3 = = 3 。x l t + 3 * x 2 e t 3 ) o u t d = ( x 1 ， t = = 4 * x l t + 2 * x 2 i t ， 9 中南工韭大学硕士学位论文 x 2 ， t = = 3 * x l t + 3 * x 2 t ) m a t h e m a t i c a 系统中的d s o l v c 函数很易给出其通解 i n 2 ：= d s o l v e r t e m p ， x l r t ，x 2 t ，t o u t r 2 = ( x l t 卜 班丝呲控里皿 丛盟里里趾， ；z - - 盥堕皿址型蔓单盟业址卿 ) 如果考虑方程( 1 2 1 ) 满足初始条件 x l ( 0 ) = 1 ，x 2 ( 0 ) = 一1( 1 2 2 ) 的解，那么可以调用m a t h e m a t i c a 系统里的函数a p p e n d 与f l a t t e n ，修改表 t e m p 的值，即： i i l 3 _ t e m p = f l a t t e n a p p e n d t e m p ， x l 0 _ = 1 x 2 o 一1 ) o u t 3 = x 1 7 t = = 4x l t + 2x 2 t ，x 2 ， t = = 3x i t + 3x 2 i t ， x 1 0 一= 1 ，x 2 0 = = 一1 ) 下面再次调用i ：x s o l v e ，这时将产生一个特解。即 i l l 4 ：= d s o l v e t e m p ，( x l t ，x 2 e t ) ，t o u t 4 一 ：。 1 3 一 隧喾 x ：i t 一 丛掣) 上面d s o l v e 稂方便地解出( 1 2 1 ) ，以及( 1 2 1 ) 与( 1 - 2 2 ) 组合问题解。另一种 方法是考虑齐次常系数系统x ，= a x ，它的通解可以写成x ：e x p ( a t ) c ，其中a 是系数矩阵，c 是一个含有初始条件的向量，对于( 1 2 1 ) 来说其系数矩阵a 为： i n 5 ：= a = “4 ，2 ) ， 3 ，3 ) o u t 5 一 】d 中南工业大学硕士学位论文 4 ，2 ) ， 3 ，3 ” 调用m a t h e m a t i c a 系统的函数m a t r i x e x p 产生矩阵幂e x p ( a t ) 。 i n 6 ：= m a t r i x e x p e a t o u t r 6 - 1 = 鲤警塑，2e r ( - - i + e s , ) ) ， 3e ( - - 。i + e s t ) , e t ( 3 - - i - 2 e s ) ) 5。 现在可以利用矩阵的点积求x = e x p ( a t ) c 。 【n 7 ：= ( 1 ，i o t e t = f - - 2e t ( - - i + e s ) 一t - 戛：1 2 ! 至兰 5 5 3e t ( 一1 十e 5 。)e ( 3 + 2e 5 ) 、 一 亏一一i 一一一7 h e b 3 ：= s i m p l i f y o u t e 8 3 = e ( 4 - f ce s ) ，e ( - - _ i 6 _ + 一e 5 ) c e j 得到的结果，除表现形武外，与第一种方法得到的结果完全一致。 1 3 摄动方法概论 在数学和工程实际问题串，经常会遇到含有小参数e 的微分方程，特别是非 线性的微分方程定解问题。小参数可能出现在微分方程中，也可能出现在定解 条件里，称这类含小参数e 的数学问题为“摄动问题”。这类问题的求解，一般采 用数值方法和近似解析方法或者两种方法兼用。按小参数e 对解进行渐近展开 是近似解法常见的方法，它是相对于小参数e = o 的退化问题的解加以校正，得 出摄动问题的近似解。 1 3 1 渐近序罗q 和渐近展开式 定义1 设 6 ( e ) ) ( n o ，1 ，2 ) 为参数e 的函数序列，如果对于每一个 1 1 中南工业大学硕士学位论文 n 0 ，当e o 时，满足条件 6 n + 1 ( e ) = o k ( 0 7 ( 1 3 1 ) 则称函数序列 6 。( ) ) 为“渐近序列”。 定义2 对于已给定函数f ( e ) 及渐近序列( 6 。( e ) ，若存在一组常数a o ，a ， ，a 。一l ，使成立下式： n 一1 f ( ) = 如k ( ) + o 6 n ( 8 ) ( 1 3 2 ) 则称n 项之和 n l za 。6 丑( ) a - d 为函数f ( o 的“n 项( 或n l 级) 渐近展开式”，其中“o ) a o ”为阶符。 定义3 若对于任意正数n ，式( 1 3 2 ) 总成立，则可记( 1 3 2 ) 式为 f ( ) a 。6 n ( 8 )( 1 3 3 ) 并称( 1 3 3 ) 式的右部为函数f ( o 的“渐近展开式”或“渐近级数”。 定义4 若f ( x ，) z 玉( x ) 6 n ( e ) ( 1 3 4 ) 对变量x 一致成立，则称( 1 3 4 ) 的右边为函数f ( x ，e ) 的“一致有效渐近展开 式”，否则称为“非一致有效渐近展开武”。 性质1 ( 1 3 3 ) 式右部的系数a 。是由f “) 唯一确定的，且有 l 一1 f ( ) 一a 丑6 n ( 8 ) 乱2 熘i 两一 性质2 若f ( x ，e ) 的渐近展开式是一致有效的，则f ( x ) k ( ) 相对于f 一- ( x ) 6 。一。( e ) 来说必须是小量。亦即厶( x ) 的奇性不比厶一z ( x ) 的奇性更强。 1 3 2 正则摄动与奇异摄动 ：在摄动问题中，根据渐近展开式的形式和e 一0 时的收敛性，将它分为正则 摄动问题与奇异摄动问题。 定义l 若摄动i 可题的解u ( x ，e ) 能够用一个e 的幂级数表示。 u ( x ，) ze 1 l 。( x ) ， n - - o 并且在区域q 内该展开武关于x 一致收敛。则称p ( x ) 为区域q 内的“正则摄 动问题”，否则称为“奇异摄动问题”。 对于正则摄动问题有三种解法：级数解法，参数解法，逐次逼近法。但用正 1 2 中南工业大学礤士学位论文 则解法，有时得到的解会失效，引起失效的主要原因有：( 1 ) 与无限域有关；( 2 ) 最高导数项有小参数 ( 3 ) 在求解过程中，微分方程的类型发生交化f ( 4 ) 出现奇 点。下面以d u f f i n g 方程为例来说明。 考虑d u f f i n g 方程 o + u + u 3 = 0( 1 3 5 ) 的一阶展开式的解 u ( t ，e ) = u o ( x ) + e u l ( x ) + 0 ( 2 )( 1 3 6 ) 将( 1 3 6 ) 代入( 1 3 5 ) 便得 程o + u o + e ( n l + u 1 + u 8 ) + 0 ( 2 ) 一0 于是便有： n o + u o = 0( 1 3 7 ) n 1 + u l + u 3 = 0( 1 3 8 ) 由( 1 3 7 ) ，( 1 3 8 ) 便得一阶近似解 u ( t ，s ) 一ac o s ( t + b ) + a 3 一毒s i n ( t + 8 ) + 蠢c o s ( 3 t + 3 b ) ( 1 3 9 ) vu - 在式( 1 3 9 ) 中，右边的第2 项出e t ，是0 ( 1 ) 阶，并不是第一项的校正项，表明正 则方法求解失效。 为了校正正则解法，到目前为止已经有多种特殊方法可以用来弥补正则解 法的失效性，如l i n d s t e d tp a i n c a r e 方法，多尺度法，平均法等，但很少见到用计 算机实现其自动求解过程。 1 4 分叉问题与化简 微分动力系统的结构稳定性讨论是研究一族相邻的动力系统的拓扑性质 之问的联系当动力系统在受到扰动后变为“邻近的”动力系统，那么在什么样 的条件下系统的拓扑结构才保持不变呢? 这就是结构稳定性问题，如果某个动 力系统是稳定的，则任意小的适当的扰动都会使系统的拓扑结构发生突然变 化，这就是分叉( 分支、分歧、分岔 。费力系统的研究不但要讨论结构稳定问题， 而且还在考虑由于结构不稳定性问题而 l 起的定性性态的可能变化。分叉问题 便成了动力系统和非线性微分方程研究的重要组成部分。 分叉问题习惯上分为静态分叉与动态分叉，静态分叉的研究不但提供平衡 点的数目随参数交化的信息，而且往往与平衡点稳定性的变化密切相关。一般 1 3 中南工业大学硕士学位论文 来说，往往是利用l i a p u n o vs c h m i d t 方法与中心流形方法，将原方程进行约化 ( 降维) 后，再等价地研究其约化方程。中心流形方法是考虑其相应的线性化算 子l 的特征值没有正实部的情形，即稳定系统开始失稳的临界状态，根据中心 流形定理( 见文献 2 ) 中心流形土的流足以反映系统的动力学性态。在本文的 最后一章用计算机代数系统实现自动化筒过程。 中南工业大学硬士学位论文 第二章奇异摄动的机器求解 在这章，应用m a t h e m a t i c a 系统本身的符号运算能力，实现了几种奇异摄动 的机器求解方法。设计了包含几个自动求解奇异摄动问题的函数的程序包；可 用它来求解一些摄动系统的有效渐近解。由于方法本身的缺陷，有时求出的解 可能是非一致有致的，可以直接从解的形式分析其有效性；有时甚至会导致错 误的解( 见文献e c 。在这里主要针对弱非线性系统求解问题的机器实现。 2 1 l i n d s t e d t _ p o i n c a r e 方法 l i n d s t e d t _ p o i n c a r e 摄动方法是获得微分扰动系统近似解析解的一种典型 的方法，是为了研究方程的近似解析解避免出现长期项，由l i n d s t e d t 等人发明 的一种方法。 2 1 1l i n d s t e d tp o i n c a r e 方法的基本思想 由第一章的( 1 3 9 ) 式知道，利用渐近级数展开法，求d u s t i n g 方程渐近展 开式解会出现长期项“t ”，从而导致渐近展开式解失效。而用l i n d s t e d tp o i n e a r e 方法可以避免长期项的出现，其基本思想是伸缩参数，使得非线性系统的频率 从线性的( 1 ) o 改变为m ( e ) 。为此引进变量交换 t = c o t 同时把。与u 按e 的幂展开为： u ；u o ( t ) + e u l ( t ) + e 2 u 2 ( t ) + = 蜘- t ：e 嘶+ 2 铆+ 下面以d u f f f n g 方程为例来说明其基本思想，考虑 a + u + e u 3 = 0 ( 2 1 1 ) 的渐近展开式解，其中n 表示对t 求两阶导数。为了防止长期项出现，设 t = 僦( 2 1 2 ) = 1 + e 蜘+ e 2 也+ ( 2 i 3 ) u = u o “) + g u l ( t ) + 8 2 1 1 2 h ) + ( 2 1 4 ) 把( 2 1 2 ) ，( 2 1 3 ) ，( 2 1 4 ) 代入( 2 1 1 ) ，并保留到的一阶为止，这样可以得 到 1 5 中南工业大学硕士学位论文 n o ”+ u o + ( u l ”+ u 1 + u 0 3 + 2 e l u 0 ”) + = 0 其中u ”表示u 对t 求两阶导数。 令e 的同阶系数相等，便得 u o + u o = 0 u l ”+ “l + 1 1 0 s + 2 0 l u o ”= 0 先求出( 2 1 6 ) 式的通解u o ，再代入( 2 1 7 ) 式求出u - 的特解，得到： u ac o s “十b ) + 专8 ( 2 ( i ) 1 a 一号a 3 ) ts i n ( t + b ) + ( 2 1 5 ) ( 2 1 6 ) ( 2 1 7 ) 蠢a 3 c o s ( 3 t 十3 p ) + ( 2 1 8 ) ( 2 1 8 ) 武的右边出现了混合长期顼s t ，会使展开武交得非一致有效，为此必须 消去混合项( 或长期项) ，这样便禄 铆= i o i l ( 2 1 9 ) q， 于是( 2 1 。8 ) 交为 u ac o s “+ p ) + 去8 a 5 c o s ( 3 x - l - 3 b ) - 1 - ( 2 1 1 d ) 是一致有效的，即满足后一项是前一项的修正。 2 1 2 l i n d s t e d t _ p o i n c a r e 方法的机器实现 这里只考虑单自由度的弱非线性系统 u + o ) 0 2 u = e f ( u ，矗) 的渐近展开式解的鸯动求解问题。 在求解过程中，首先引进变量变换( 2 1 2 ) 式，参数。与u 关于小参数e 的 幂级数形式，即( 2 1 3 ) 式与( 2 1 4 ) 式。应用m a t h e m a t i c a 的替换规则，将 ( 2 1 2 ) ，( 2 1 3 ) ，( 2 1 4 ) 式代入原方程，展开所有项并令8 同次幂项相等，这 样可以得到一系列的微分方程。再迭代求解，首先解出u o ，然后把u o 代入到其他 的方程中，并把正弦函数与余余弦函数的幂次项或两者乘积项化简为正弦函数 或索弦函数的倍角形式，这样便于消去长期项。如此进行，可以求得任意阶的展 开武解。应用m a t h e m a t l c a 系统中的替换规则，匹配功能，s o l v e 函数，d s o l v e 函数可以在微机上实现自动求解，其运行过程可以大体描述如下t l 。输入微分方程e x p r i - - - - - e x p r 2 ，小参数e p s u o n ，未知函数符号u e t j ，频率 f _ ) 。，时间变量t ，| 羚数n ， 1b 中南工业大学硕士学位论文 2 。作变量替换t a u c o t ，令= 1 - e p s i l o n c 0 1 n u e p s f l o n 2 - ( i ) 2 + ( 若原方 程不显含小参数e p s f l o n ，则c o = e p s i o n t o l - - e p s f l o n 2 2 + ) ，u u o t + e p s f l o n u li t - - e p s i l o n 2 u 2 t + ，化简微分方程； ， 3 。按e p s i l o n 的幂合并同类项，取出e p s f l o n 的同次幂组成n + 1 个微分方 程 4 。求出u d 的通解，当k l 时，求u 。的特解 5 。把求出的u ”u ：，u t 代入下一个方程，进行三角函数化简，消去长期 项，解出啦+ ，。k = k + 1 ，并判断是否超过指定的阶数，若未超过转到步骤4 。， 酽输出u 的渐近展开式与m 的表达式。 下面是关于主函数和两个子函数的说明。 ( 1 ) l i n d s t e d t 函数的输入形式为t l i n d s t e d t e x p r l = = e x p r 2 ， e p s f l o n ，n ，( u t ，t ，w o ) 。 输出形式为： “o m e g a - * ) ， u t 肌 一) ) 其中o m e g a 为频率参数，输出结果是一个表，表的第一项是参数交换的表达式， 第二项为所求渐近近似解的表达式。 ( 2 ) t r i $ s i m p l i 印 数是化筒表达式中的三角函数，园为每次迭代的时 候，都需要分离长期项与短期项，这样必须把三角函数s i n x ，c o s b 的乘积形 式或乘幂形式变成积化和差的形式( 倍角的形式) ，即如下列那样： ( s i n t 弘；寺( 1 一c 销 2 t ) ： ( c o s t ) 2 一告( 1 + c 2 t ) ， 11 s i n t 1 * c o s t 2 = 告s m t l t 2 一号s h t 1 + t 2 等。 ( 3 ) s p e c i a l i ) s o l v e 薮用来求形如l u 一 t + c o 口2 u t = f c a s a t ，s i n b t 的特解，如u 一 t + u t = c 3 t ，利用s p e c i a l d s o l v e 求解后的结果为： u t 一 一告c o s 3 t ) ) 。 2 1 3 应用攀例与分的 这里考查两个例子，检查一下l i n d s t e d t 函数的运行情况，首先看一看 1 7 中南工业大学硬士学位论文 d u f f i n g 方程 n + u + u 3 = 0 的有效渐近展开式解。 i n i ：2 十u c l l l i n d s t e d tu ” t + e * u - t “3 = 0 ，( e ，1 ) ， u t ，t ，i ) o u t 1 = t a u - - o m e g at , o m e g a - - 1 + 羔皇雩警， u t 一 a l p h ac o s b e t a - b t a u + 型坚旦嘎譬划) 在p e n t i u ms c p ui o o m h z ，内存为i b m 的微机上的运行时间约为6 秒， 结果与文献 4 3 相符。 函数l i n d s t e d t 可以计算出精确到任意阶的有效渐近近似解。但是随着所 要精确的阶数增加，出现的未知函数就增加，每一步迭代所要处理的已知函数 的微分，三角函数的化简会急剧增多，亦即m a t h c m a t i c a 系统所调用d 函数， d e r i v e t i v e 函数，调用子程序t r i g s i m p l i f y 会急剧增加，而且调用的系统函 数d s o l v e 所处理的表达式会随之交得很复杂，这样，运算的复杂性就很大，例 如上面的d u f f i n g 方程的渐近解，当要求精确到三阶时。在p e n t i u msc p u i o o m h z ，内存为1 6 m 的微机上的运行时问约为1 分2 6 秒。作为改进方案，主要 避开调用系统d s o l v e 口函数求第二步迭代的每个方程的特解，改进方案见第四 章。 第二个例子，考虑微分方程： 诅+ 4 u e 亡2 u = oe 1 + 型旦；址， u t 一 a l p h ac o s b e t a 十2t a u 一a l p h a a e c o s 、3 ，b e t a + 6t a u ) 在p e n t i u m sc p u1 0 0 m h z ，内存为1 6 m 的微机上运行时间约为4 秒，作为 检验其结果的正确性，调用m a t h e m a t i c a 系统中n d s o l v e 函数求( 2 i 1 1 ) 式 1 8 中南工业大学硬士学位论文 的数值解，并调用p l o t 函数画出数值解的图形与用l i n d s t c d t 函数所求解 ( 2 1 1 2 ) 式的图形，见下面的图2 1 1 ( 一1 ) ，图2 1 2 ( e = 0 5 ) ，图2 1 3 ( 一 0 1 ) 。图示表明当e 较小时，用l i n d s t e d t 函数与用数值解法的结果非常接近。 八瓜一i 1 v4 蝴 一0 一o 0 0 一 o 一 0 图2 1 1 摄i 由衅 教膻畔 八一“ l j 图2 1 2 。、厂、一 卜“l 1 v 4v 图2 1 3 最后一个例子考虑v a nd e rp o l 方程 n + u + e ( u 2 - 1 ) f i ；0 的一阶渐近展开式解 1 9 摄却筒耳 数值衅 中南工业大学硕士学位论文 i n 4 ；一 l i n d s t e d t u ” t + u e t + o ( u f t 7 2 一1 ) * u t t = = 0 ，( e ，1 ) ，( u e t ， t ，1 ) o u t 4 = t a u 一 o m e g a * t ，钿e g a 一 1 + e * o m e g a l ， u t 一 a l p h a * c o s b e t a + t a u + ( a l p h a * e * ( 4 * i * c o s b e t a + t a u 一 i , a l p h a 2 * c o s b e t a + t a u + 1 6 * t a u * c o s b e t a + t a u 一 4 * a l p h a 2 * t a u * c o s b e t a + t a u 一4 * s i n b e t a + t a u + 一一a l p h a 2 * s i n b e t a ! t a u - a l p h a ：2 * s _ i n 3 * ( b e t a 一+ t a _ 幽! ) 3 啦 ( 2 1 1 3 ) 在p e n t i u m sc p u1 0 0 m h z ，内存1 6 m 的微机上运行时间为2 9 秒，结果与 文献e s 相符。 但是结果( 2