求解Bessel方程.doc

1、1一类 Bessel 方程在不同边界条件下的差分解信息与计算机科学专业 赵鑫指导老师：李顺初 教授【摘要】 本文主要 阐述的是使用差分法解一类 Bessel 方程在不同边界条件下的离散解，以及使用追赶法求解得到的线性方程组.在得到方程的离散解后，利用 MATLAB 程序，绘制方程的离散解，并 分析这些离散解构成的解曲线，及各种参数的变化 对解曲线的影响。关键字：微分方程 Bessel 方程 差分法 MATLAB 程序 The kind of Bessel equation dissolve in different boundary【Abstract】 This thesis introduc

2、es solving the kind of Bessel equation in different boundary by difference method and gains discrete solution, and solving the linear equation by chase method . Protracting the discrete solution of the Bessel equation by Matlab, and observing the curve which is compose of the discrete solution, disc

3、ussing the condition of the curve in different parameter. Key word: Bessel function Matlab program Difference method一 前言Bessel 函数是工程技术中常用的一种特殊函数，主要地来源于 圆柱形或圆柱物体的有关物理问题，因此也叫做圆柱函数。这种函数是研究 一些物理问题时所归结 成的贝塞尔方程的解。022 xyxyx早在贝塞尔之前，欧拉（Euler）、拉格朗日（ Lagrange）等人研究过这个方程，并求出了它的通解。贝塞尔在1824 年关于天文学问题的研究中，系 统地研究了此方程

4、的解，同时编制出这些函数具有十位小数的函数表。此后，称这种函数为贝塞尔函数，并被广泛地应用到物理学和技 术科学中。本文主要阐述的是使用差分法解一类 Bessel 方程在不同边界条件下的定解问题的离散解，并借助MATLAB 程序的绘图工具对 其给予直观的表现，分析在不同边界条件下参数的变化对解曲线的影响。二 函数模型及其求解2.1 模型的建立本文要求解的模型如下： .1,0,0)(0112 RuRbayoruQdxybax 均 为 常 数 ； 且其 中 ： 右 边 界 条 件左 边 界 条 件22.2 模型的分析与求解采用差分法求解上面的模型，首先应该利用导数的定义将 导数化为：(1)hxffh

5、xffdxyh )(lim0 (2)2002 )()( )()()( ) lilihxffxf hxffff hh 先写出差分格式。将（1）（2）代入微分方程 得到：012uydxy(3) )2()()()( hfhfxhfxuh在区间1，R上，将这个区间等分成 n 等份，每一等份的大小为 ，为了求出解曲线上nR/)1(的这 个点的 值，故 分别 取为： .代入(3) 式得到一个含 个未知数线性方1nxn1,.21, 程组: 0)1()2(1)(1)2(1)( .323)()( . 0)51()3()41(5)31(1)( 42)32(2( 0)()(22 nhfnhnfhnhfu fff h

6、fhfhfhu ffff.(4)方程组(4)含有 个方程,为了解出这个线性方程组,必须再增加两个条件,即左边界和右边界条件,这里先列出左边界条件的情况,将(1)式代入)(|1uqdxyba得到: hhffyx1|)(即: (5)qbffbah)()1(3右边界有两种情况,下面将分别讨论.2.2.1 第一种右边界条件 0)(Ry因为 所以nh1(6)(故由(4)(5)(6)可以得到一个 阶线性方程组:n 0)1()2(1)2(1)2(1)( .03)32. )51()3()41()52()31()1( 04220)()()(1(2 hnfhnhnfhnuh fff hfhfhfhu ffffqh

7、bffah(7)2.2.2 第二种右边界条件 0)(Ry根据导数的定义得: 0)(|)(| hRffxfhfdxyRxR即: )(ff也即: (8)11nhfn故由(4)(5)(8)可以得到一个 n 阶线性方程组: 0)1()2(1)2(1)2(1)( .0)32 )3()()( . )51()3()41()52()1()1( 04232 )()()()()()( 02111)()(2 hnfhnhnfhnuh fff hfhfhfhu fff hhfhuqbffa4(9)在 Matlab 程序实现中将采用追赶法来解出 即为微分方程解).1().21(),1( hnfhff 曲线上的点.三 M

8、atlab 程序实现31 Matlab 工具 简介MATLAB 是一种功能非常强 大的科学计算软件。它源于 MATrix LABoratory 一词，愿意为矩阵实验室。一开始它是一种专门用于矩阵数值计算的软件。随着 MATLAB 逐渐市场化，它不仅具有了数值计算功能，而且具有了数据可视化功能。在 MATLAB6.x 中，它不 仅在数 值计算、符号运算和 图形处理功能上进一步加强，而且又增加了许多工具箱。MATLAB 可以说自产生之日起，就以其强大的功能和良好的开放性而在科学计算诸软件中独占熬头。如今，新版本的 MATLAB 在数值计算、符号运算及图形处理方面都在同类产品中占有优势。MATLAB

9、 4.0 以上（不包括 4.0 版本）的各版本，不 仅在数值计算上继续保持着相对其他同类软件的绝对优势，而且还开发了自己的符号运算功能。特别是 MATLAB6.x 版本在符号运算功能上丝毫不逊色于其他各类软件，如 MathCAD, Mathematica 等。 这样，用 户就不用像以前的 计算机人员那样在掌握 MATLAB 的同时还有学习另一种符号运算软件。用 户只要学会了 MATLAB6.x，就可以方便地处理诸如矩阵变换及运算、多 项式运算、微积分运算、线性与非线性方程求解、常微分方程求解、偏微分方程求解、插值与拟合、统计及优化等问题了。MATLAB 还允许用户以数学形式的 语言编写程序，比

10、 BASIC、C 语言和 FORTRAN 等语言更接近于书写计算公式的思维方式。它的操作和功能函数指令就是以平 时计算机和数学书上的一些简单英文单词表达的。由于它在很长一段时间内是用 C 语言开发的，它的不多的几个程序流控制语句同 C 语言差别甚微，初学者很容易掌握。从形式上看，MATLAB 程序文件是一个纯文本文件， 扩展名为 m。用任何字处理软件都可以对它进行编写和修改，因此程序易调试，人机交互性强。本文借助 MATLAB 强大的绘图 工具，直观的描述了在参数变化时函数曲线的变化。3.2 算法及算法分析3.2.1 追赶法在用差分法求解二阶微分方程的边值问题时，常常都会遇到三对角线形方程组

11、AX=f,即：=nnbaccba11221. xx21.nff121.其中系数矩阵 A 的元素满足条件：(10)niiiiabcc001且5这类方程组是适合于用三角分解法求解的典型问题之一，根据系数矩阵 A 的特点，进行 LU 分解：A= =LUnn121. 1.1其中 为待定系数，比 较矩阵 A 与 LU 的对应元素，有：ii,.(11)1.2.3,1111nicbiaciiii由关系式（10）（11）可以看出： 011cb故 有意义，且 ；进而又可以看出：1/c0221212122 cababa故 有意义且 ，可用 归纳法证明：/c0（即： ） i=1，2，n-1iiicbi因此可以从(1

12、1)式解出待定数niabciii iiiii .3,2/.111另一方面,由 及 知 不等于 0,故矩阵 L 非奇异,从而满足条件(11)的三对角线0i na矩阵 A=LU 非奇异,方程组 AX= 有唯一解.f用上述 LU 分解法解线性方程 组的过程可以归纳为:(1) 实现 A=LU 的分解,按递推式计算 .ii,(2) 求解方程组 LY= ,相应的递推算式是: b niabyafyiiiii .3,2)/()(/ 111(3) 求解方程组 UX=Y,相应的递推算式是:6nixyxiiin .3,21由于计算 及 的过程为追的过程，计算方程组的解21.n nyy.1的过 程称为赶的过程，故上述

13、求解方程组的方法称为追赶法。1.xxn下面将采用追赶法分别讨论两种不同右边界条件下线形方程组的解，即线形方程组(4)与(9)的解。3.2.2 第一种右边界条件 0)(Ry（7）式的系数矩阵为: )2(1)2(1)0.00 33)3(1. 0012)1( 0.)1( 222 hnhnuhnuhhuhba若要使方程组有唯一解,则方程组的系数矩阵必须满足下列条件:(12) 0 3.0)1)(1)(2 202u nihiuhiihnnbah条件(12)的作用是使线性方程组(7)的系数矩阵非奇异,从而保证线性方程组有唯一解,再按照追赶法的(1)(2)(3)步求解线性方程组。解方程组(7)的追赶法算法的流

14、程图见附录 1.程序的源代码见附录 2:3.2.3 第 2 种右边界条件 0)(Ry（9）式的系数矩阵为: )2(1)2(1)0.00 33)3(1. 0012)1( 0.)1( 222 hnhnuhnuhhuhba若要使方程组有唯一解,则方程组的系数矩阵必须满足下列条件:7 0 3.0)1(1)(22)02u nihiuhiihnnbah解方程组(7)与(9) 的算法相同,流程图也相同.解方程组(9)的源代码见附录 3:四 模型的分析与结果根据前面的求解，我们得到了 这个模型的离散解。下面将 讨论 在不同边界条件下,不同参数 a,b,u,q 对解曲线的影响,并绘制图版,其中横坐标是 x 轴，纵坐标是 y 轴.4.1 参数 u 对解曲线的影响4.1.1 第一种边界条件 0)(Ry84.1.2 第 2 种边界条件 0Ry9104.2 参数 a 对解曲线的影响4.2.1 第一种边界条件 0)(Ry