偏微分方程数值解法
- 格式:doc
- 大小:415.50 KB
- 文档页数:6
偏微分方程的数值方法偏微分方程(Partial Differential Equations,简称PDEs)是数学中研究的重要分支,广泛应用于物理学、工程学等领域中。
由于一些复杂的PDEs难以找到解析解,因此需要借助数值方法进行求解。
本文将介绍偏微分方程的数值解法,包括有限差分法、有限元法和谱方法等。
一、有限差分法(Finite Difference Method)有限差分法是解偏微分方程最常用的数值方法之一。
它将偏微分方程中的导数用差商来近似,将空间离散成若干个小区间和时间离散成若干个小时间步长。
通过求解离散化后的代数方程,可以得到原偏微分方程的数值解。
以二维的泊松方程为例,偏微分方程可以表示为:∂²u/∂x² + ∂²u/∂y² = f(x, y)其中,u(x, y)为未知函数,f(x, y)为已知函数。
我们可以将空间离散成Nx × Ny个小区间,时间离散成Nt个小时间步长。
利用中心差分法可以近似表示导数,我们可以得到离散化的代数方程组。
二、有限元法(Finite Element Method)有限元法是一种重要的数值解PDEs的方法。
它将求解区域离散化成一系列的单元,再通过插值函数将每个单元上的未知函数近似表达。
然后,利用加权残差方法,将PDEs转化成代数方程组。
在有限元法中,采用形函数来近似未知函数。
将偏微分方程转化为弱形式,通过选取适当的形函数和权函数,可以得到离散化后的代数方程组。
有限元法适用于求解各种各样的偏微分方程,包括静态和动态、线性和非线性、自由边界和固定边界等问题。
三、谱方法(Spectral Method)谱方法是一种基于特殊函数(如正交多项式)的数值方法,用于解PDEs。
谱方法在求解偏微分方程时,利用高阶连续函数拟合初始条件和边界条件,通过调整特殊函数的系数来近似求解解析解。
谱方法具有高精度和快速收敛的特点,适用于各种偏微分方程求解。
偏微分方程数值解法
偏微分方程数值解法是一种利用计算机技术获取偏微分方程数值解的方法,它主要目标是解决微分方程的精确、快速、可靠的数值解。
偏微分方程数值解法交叉应用于分析数学、力学、电磁学等不同领域的各种模型,能够大大提高解决微分方程的效率。
偏微分方程数值解法大致分为两个方面:一是求解偏微分方程的离散数值解法;二是精确解对分解数值解法,如多阶谱方法、牛顿法和共轭梯度法等。
其中,离散数值解法是把偏微分方程抽象成一系列数值求解问题,并进行递推叠加求解,而精确解对分解数值解法则是通过优化问题方式求解微分方程精确解,以达到精确求解的目的。
偏微分方程数值解法的有效解决的方法,给科学与技术研究带来了很大的帮助。
它不但克服了无法精确解决某些复杂偏微分方程的困难,而且有更快的求解效率,也可以很好地满足实际科技应用的需要。
偏微分方程数值解法的应用已经普遍发挥出重要的作用,不仅可以解决物理科学问题,还可以解决经济学、商业投资、财务分析等复杂的数学模型。
因此,偏微分方程数值解法的应用已在各个领域得到了广泛的应用,为科学与技术研究提供了很大的帮助,在微分方程求解方面产生了重要的影响。
第十六章 偏微分方程的数值解法科学研究和工程技术中的许多问题可建立偏微分方程的数学模型。
包含多个自变量的微分方程称为偏微分方程(partial differential equation),简称PDE 。
偏微分方程问题,其求解是十分困难的。
除少数特殊情况外,绝大多数情况均难以求出精确解。
因此,近似解法就显得更为重要。
本章仅介绍求解各类典型偏微分方程定解问题的差分方法。
16.1 几类偏微分方程的定解问题一个偏微分方程的表示通常如下:(,,,,)x x x y y y x y A B C f x y Φ+Φ+Φ=ΦΦΦ (16.1.1) 式中,,,A B C 是常数,称为拟线性(quasilinear)数。
通常,存在3种拟线性方程: 双曲型(hyperbolic)方程:240B AC ->; 抛物线型(parabolic)方程:240B AC -=; 椭圆型(ellliptic)方程:240B AC -<。
16.1.2 双曲型方程最简单形式为一阶双曲型方程:0u ua t x∂∂+=∂∂ (16.1.2) 物理中常见的一维振动与波动问题可用二阶波动方程:22222u u a t x∂∂=∂∂ (16.1.3) 描述,它是双曲型方程的典型形式。
方程的初值问题为:2222200,(,0)()()t u uat x tx u x x u x x t ϕψ=⎧∂∂=>-∞<<+∞⎪∂∂⎪⎪=⎨⎪∂⎪=-∞<<+∞⎪∂⎩ (16.1.4)边界条件一般有三类,最简单的初边值问题为:2222212000,0(,0)(0,)(),(,)()0()t u ua t T x l t x u x lu t g t u l t g t t T ux x t ϕψ=⎧∂∂==<<<<⎪∂∂⎪⎪=≤⎪⎨==≤≤⎪⎪∂=-∞<<+∞⎪∂⎪⎩ (16.1.5)16.1.3 抛物型方程其最简单的形式为一维热传导方程:220(0)u ua a t x∂∂-=>∂∂ (16.1.8) 方程可以有两种不同类型的定解问题:(1) 初值问题:2200,(,0)()u ua t x t xu x x x ϕ⎧∂∂-=>-∞<<+∞⎪∂∂⎨⎪=-∞<<+∞⎩(16.1.6)(2) 初边值问题:221200,0(,0)()0(0,)(),(,)()0u ua t T x l t x u x x x l u t g t u l t g t t Tϕ⎧∂∂-=<<<<⎪∂∂⎪⎪=≤≤⎨⎪==≤≤⎪⎪⎩(16.1.7) 其中()x ϕ,1()g t ,2()g t 为已知函数,且满足连接条件:12(0)(0),()(0)g l g ϕϕ== (16.1.8)边界条件12(0,)(),(,)()u t g t u l t g t ==为第一类边界条件。
偏微分方程组数值解法
偏微分方程组是描述自然、科学和工程问题的重要数学工具。
由于解析解通常难以获得,因此需要使用数值方法来解决这些方程组。
本文将介绍偏微分方程组的一些数值解法,包括有限差分法、有限元法、谱方法和边界元法等。
有限差分法是一种基本的数值方法,将偏微分方程转化为差分方程,然后使用迭代算法求解。
该方法易于理解和实现,但对网格的选择和精度的控制要求较高。
有限元法是目前广泛使用的数值方法之一,它将偏微分方程转化为变分问题,并通过对函数空间的逼近来求解。
该方法对复杂几何形状和非线性问题有很好的适应性,但需要对网格进行精细的划分,计算量较大。
谱方法是一种高精度的数值方法,它将偏微分方程转化为特征值问题,并使用级数逼近来求解。
该方法在高精度求解、解析性质研究和数值计算效率方面具有优势,但需要对函数的光滑性和周期性有较高的要求。
边界元法是一种基于边界积分方程的数值方法,它将偏微分方程转化为边界积分方程,并使用离散化方法求解。
该方法适用于求解边界问题和无穷域问题,但对边界的光滑性和边界积分算子的性质有较高的要求。
总之,在实际问题中选择合适的数值方法需要综合考虑问题的性质、计算资源、精度要求等因素。
偏微分方程数值解法偏微分方程(Partial Differential Equations,简称PDE)是数学中重要的研究对象,其在物理学、工程学、经济学等领域有广泛的应用。
然而,对于大多数偏微分方程而言,很难通过解析方法得到精确解,因此需要借助数值解法来求解。
本文将介绍几种常见的偏微分方程数值解法。
一、有限差分法(Finite Difference Method)有限差分法是一种常见且直观的偏微分方程数值解法。
其基本思想是将偏微分方程中的导数通过差分近似来表示,然后通过离散化的方式转化为代数方程组进行求解。
对于一维偏微分方程,可以通过将空间坐标离散化成一系列有限的格点,并使用中心差分格式来近似原方程中的导数项。
然后,将时间坐标离散化,利用差分格式逐步计算每个时间步的解。
最后,通过迭代计算所有时间步,可以得到整个时间域上的解。
对于二维或高维的偏微分方程,可以将空间坐标进行多重离散化,利用多维的中心差分格式进行近似,然后通过迭代计算得到整个空间域上的解。
二、有限元法(Finite Element Method)有限元法是另一种重要的偏微分方程数值解法。
其基本思想是将求解区域分割成有限数量的子区域(单元),然后通过求解子区域上的局部问题来逼近整个求解区域上的解。
在有限元法中,首先选择适当的形状函数,在每个单元上构建近似函数空间。
然后,通过构建变分问题,将原偏微分方程转化为一系列代数方程。
最后,通过求解这些代数方程,可以得到整个求解区域上的解。
有限元法适用于各种复杂的边界条件和几何构型,因此在实际工程问题中被广泛应用。
三、谱方法(Spectral Methods)谱方法是一种基于特定基函数(如切比雪夫多项式、勒让德多项式等)展开解的偏微分方程数值解法。
与有限差分法和有限元法不同,谱方法在整个求解区域上都具有高精度和快速收敛的特性。
在谱方法中,通过选择适当的基函数,并利用其正交性质,可以将解在整个求解区域上展开为基函数系数的线性组合。
偏微分方程是数学中的一大重要分支,广泛应用于物理、工程、金融等领域。
其求解方法可以分为解析解法和数值解法。
解析解法要求方程具有可积性,适用于一些简单的方程,但是对于复杂的方程往往无法得到解析解。
而数值解法通过将方程离散化,利用数值计算方法得到数值解,是一种弥补解析解法不足的重要手段。
在高等数学中,偏微分方程数值解法主要包括差分法、有限元法和有限差分法。
其中,差分法是最早应用于求解偏微分方程的数值方法之一。
差分法通过将偏微分方程中的导数用差商的形式来近似表示,将连续的问题转化为离散的问题,再通过计算机程序来进行求解。
差分法的优点是简单易懂、计算速度快,适用于一些较为简单的偏微分方程。
但是差分法的精度受到离散化步长的影响,不适用于一些对精度要求较高的问题。
有限元法是一种更为广泛应用的偏微分方程数值解法。
有限元法通过将求解区域分割成有限多个小区域,用简单形状的基函数来逼近真实解,再通过求解线性方程组得到数值解。
有限元法的优点在于适用于复杂的几何形状、能够处理不规则的边界条件,并且精度较高。
有限元法还具有较好的可扩展性,可以处理大规模的求解问题。
因此,有限元法在工程领域的应用非常广泛。
有限差分法是一种通过计算导数来逼近微分方程的数值解法。
有限差分法基于泰勒展开公式,将微分算子在某点处的展开为有限多个导数的差商的线性组合。
通过将微分算子离散化,可以将偏微分方程转化为代数方程组,再通过求解方程组来得到数值解。
有限差分法的优点在于简单易懂,计算速度较快。
但是由于差商的导数逼近误差,有限差分法的精度受到离散化步长的影响,需要选择合适的步长来保证精度。
总的来说,高等数学中的偏微分方程数值解法是研究偏微分方程数值计算的一大热点和难点。
不同的数值方法适用于不同的问题,需要根据具体情况来选择适合的数值方法。
在求解偏微分方程时,还需要注意数值误差对结果的影响,并通过适当选择离散化步长和网格数量等参数来提高数值解的精度。
随着计算机技术的发展,偏微分方程数值解法将会越来越广泛地应用于实际问题的求解中。
偏微分方程的数值解法偏微分方程(Partial Differential Equation,PDE)是描述物理、化学、工程学等许多科学领域中变化的方程。
由于PDE的求解通常是困难的,因此需要使用数值方法。
本文将介绍偏微分方程的数值解法。
一般来说,求解PDE需要求得其解析解。
然而,对于复杂的PDE,往往不存在解析解,因此需要使用数值解法求解。
数值解法可以分为两类:有限差分法和有限元法。
有限差分法是将计算区域分成网格,利用差分公式将PDE转化为离散方程组,然后使用解线性方程组的方法求解。
有限元法则是将计算区域分成有限数量的单元,每个单元内使用多项式函数逼近PDE的解,在单元之间匹配边界条件,得到整个区域上的逼近解。
首先讨论有限差分法。
常见的差分公式包括前向差分、后向差分、中心差分等。
以一维热传导方程为例,其偏微分方程形式为:$$ \frac{\partial u}{\partial t}=\frac{\partial^2 u}{\partial x^2} $$其中,$u(x,t)$表示物理量在时刻$t$和位置$x$处的值。
将其离散化,可得到:$$ \frac{u(x_i,t_{j+1})-u(x_i,t_j)}{\Delta t}=\frac{u(x_{i+1},t_j)-2u(x_i,t_j)+u(x_{i-1},t_j)}{\Delta x^2} $$其中,$x_i=i\Delta x$,$t_j=j\Delta t$,$\Delta x$和$\Delta t$分别表示$x$和$t$上的网格大小。
该差分方程可以通过简单的代数操作化为:$$ u_{i,j+1}=u_{i,j}+\frac{\Delta t}{\Delta x^2}(u_{i+1,j}-2u_{i,j}+u_{i-1,j}) $$其中,$u_{i,j}$表示在网格点$(x_i,t_j)$处的数值解。
由于差分方程中一阶导数的差分公式只具有一阶精度,因此需要使用两个网格点来逼近一阶导数。
§4 偏微分方程的数值解法一、 差分法差分法是常用的一种数值解法.它是在微分方程中用差商代替偏导数,得到相应的差分方程,通过解差分方程得到微分方程解的近似值. 1. 网格与差商在平面 (x ,y )上的一以S 为边界的有界区域D 上考虑定解问题.为了用差分法求解,分别作平行于x 轴和y 轴的直线族.⎩⎨⎧====jhy y ihx x i i (i ,j =0,±1,±2,…,±n ) 作成一个正方形网格,这里h 为事先指定的正数,称为步长;网格的交点称为节点,简记为(i ,j ).取一些与边界S 接近的网格节点,用它们连成折线S h ,S h 所围成的区域记作D h .称D h 内的节点为内节点,位于S h 上的节点称为边界节点(图14.7).下面都在网格D h + S h 上考虑问题:寻求各个节点上解的近似值.在边界节点上取与它最接近的边界点上的边值作为解的近似值,而在内节点上,用以下的差商代替偏导数:()()[]()()[]()()()[]()()()[]()()()[]y x u h y x u y h x u h y x u hy x u h y x u y x u h y x u hy u y h x u y x u y h x u h x u y x u h y x u hyu y x u y h x u h x u ,),(,,1,,2,1,,2,1,,1,,122222222++-+-+≈∂∂∂-+-+≈∂∂-+-+≈∂∂-+≈∂∂-+≈∂∂注意, 1︒ 式中的差商()()[]y x u y h x u h ,,1-+称为向后差商,而()()[]y h x u y x u h,,1--称为向前差商,()()[]y h x u y h x u h,,21--+称为中心差商.也可用向前差商或中心差商代替一阶偏导数.2︒ x 轴与y 轴也可分别采用不同的步长h ,l ,即用直线族⎩⎨⎧====jhy y ihx x j i (i,j =0, ±1, ±2 , ) 作一个矩形网格.图14.72. 椭圆型方程的差分方法[五点格式] 考虑拉普拉斯方程的第一边值问题()()⎪⎪⎩⎪⎪⎨⎧=∈=∂∂+∂∂y x u D y x y ux u S ,,02222μ式中μ(x ,y )为定义在D 的边界S 上的已知函数.采用正方形网格,记u (x i ,y j )=u ij ,在节点(i ,j )上分别用差商 u u u h u u u h i j ij i j i j ij i j -+-+-+-+11211222,,,,,代替2222,yux u ∂∂∂∂,对应的差分方程为u u u h u u u hi j ij i j i j ij i j -+-+-++-+=112112220,,,, (1) 或u u u u u ij i j i j i j i j =+++-+-+141111,,,,即任一节点(i ,j )上u ij 的值等于周围相邻节点上解的值的算术平均,这种形式的差分方程称为五点格式,在边界节点上取()()()h j i ij S j i y x u ∈=,,**μ(2) 式中(x i *,y j *)是与节点(i ,j )最接近的S 上的点.于是得到了以所有内节点上的u ij 值为未知量的若干个线性代数方程,由于每一个节点都可列出一个方程,所以未知量的个数与方程的个数都等于节点的总数,于是,可用通常的方法(如高斯消去法)解此线性代数方程组,但当步长不很大时,用高斯消去法将会遇到很大困难,可用下面介绍的其他方法求解. 若h →0时,差分方程的解收敛于微分方程的解,则称差分方程为收敛的.在计算过程中,由于进行四则运算引起舍入误差,每一步计算的舍入误差都会影响以后的计算结果,如果这种影响所产生的计算偏差可以控制,而不至于随着计算次数的增加而无限增大,则称差分方程是稳定的.[迭代法解差分方程] 在五点格式的差分方程中,任意取一组初值{u ij },只要求它们在边界节点(i ,j )上取以已知值μ(x i *,y j *),然后用逐次逼近法(也称迭代法)解五点格式:()()()()()[]() ,2,1,0411,1,,1,11=+++=+-+-+n u u u u u n j i n j i n j i n j i n ij 逐次求出{u ij (n )}.当(i+1,j ),(i -1,j ),(i ,j -1),(i ,j+1)中有一点是边界节点时,每次迭代时,都要在这一点上取最接近的边界点的值.当n →∞时,u ij (n )收敛于差分方程的解,因此n 充分大时,{u ij (n )}可作差分方程的近似解,迭代次数越多,近似解越接近差分方程的解.[用调节余数法求节点上解的近似值] 以差商代替Δu 时,用节点(i+1,j ),(i -1,j ),(i ,j+1),(i ,j -1)上u 的近似值来表示u 在节点(i ,j )的值将产生的误差,称此误差为余数R ij ,即()()()()()ij j i j i j i j i j i R y x u h y x u h y x u y h x u y h x u =--+++-++,4,,,,设在(i ,j )上给u ij 以改变量δu ij ,从上式可见R ij 将减少4δu ij ,而其余含有u (x i ,y j )的差分方程中的余数将增加δu ij ,多次调整δu ij 的值就可将余数调整到许可的有效数字的范围内,这样可获得各节点上u (x ,y )的近似值.这种方法比较简单,特别在对称区域中计算更简捷.例 求Δu =0在内节点A ,B ,C ,D 上解的近似值.设在边界节点1,2,3,4上分别取值为1,2,3,4(图14.8) 解 记u (A )=u A ,点A ,B ,C ,D 的余数分别为-4u A + u B + u c +5=R A u A -4 u B + u D +7=R Bu A-4 u c + u D +3=R Cu B + u c -4u D +5=R D以边界节点的边值的算术平均值作为初次近似值,即u A (0)=u B (0)=u C (0)=u D (0)=2.5则相应的余数为:R A =0, R B =2, R C = -2, R D =0最大余数为±2.先用δu C =-0.5把R C 缩减为零,u C 相应地变为2,这时R A , R D 也同时缩减(-0.5),新余数是R A =-0.5,R B =2,0=C R , R D =-0.5.类似地再变更δu B =0.5,从而 u B 变为3,则得新余数为0====D C B A R R R R .这样便可消去各节点的余数,于是u 在各节点的近似值为:u A =2.5, u B =3, u C =2, u D =2.5现将各次近似值及余数列表如下:次数调 整 值第n 次近似值及余数u A R A u B R B u C R C u D R D 0 1 2δu C = -0.5 δu B = 0.5 2.5 2.5 2.5 0 -0.5 0 2.5 2.5 3 2 2 0 2.5 2 2 -2 0 0 2.5 2.5 2.5 0 -0.5 0 结果近似值2.5322.5[解重调和方程的差分方法] 在矩形D (x 0≤x ≤x 0+a ,y 0≤y ≤y 0+a )中考虑重调和方程024*******=∂∂+∂∂∂+∂∂=yuy x u x u u ∆取步长h an=,引直线族图14.8⎩⎨⎧+=+=jh y y ihx x 00 (i , j = 0, 1, 2,, n ) 作成一个正方形网格.用差商代替偏导数()()()()()[]{()()()()[]()()()()[]}h y x u h y x u y h x u y h x u h y h x u h y h x u h y h x u h y h x u h y x u h y x u y h x u y h x u y x u 2,2,,2,2,,,,2,,,,8201,-+++-++---++-+-++++--+++-++=上式表明了以(x ,y )为中心时,u (x ,y )的函数值与周围各点函数值的关系,但对于邻近边界节点的点(x ,y ),如图14.9中的A ,就不能直接使用上式,此时将划分网格的直线族延伸,在延伸线上定出与边界距离为h 的点,称这些点为外邻边界节点,如图14.9以A 为中心时,点E ,C 为边界节点,点J ,K 为E ,C 的外邻边界节点,用下法补充定义外邻边界节点J 处函数的近似值u J ,便可应用上面的公式.1︒ 边界条件为()()()S P P x uP u SS ∈==21,μ∂∂μ 时,定义u J =u A -2μ2(E )h .2︒ 边界条件为()()()S P P xuP u SS ∈=∂∂=2221,μμ时,定义u J =2μ1(E )-u A -h 2μ2(E ). [其他与Δu 有关的网格] 1︒ 三角网格(图14.10(a ))取P 0(x ,y )为中心,它的周围6个邻近节点分别为:()()⎪⎪⎭⎫⎝⎛-+⎪⎪⎭⎫ ⎝⎛---⎪⎪⎭⎫⎝⎛+-⎪⎪⎭⎫⎝⎛+++h y h x P h y h x P y h x P h y h x P h y h x P y h x P 23,2,23,2,,23,223,2,,654321则 R u h u u u h i i +∆+∆=⎪⎭⎫⎝⎛-∑=226102161632式中u i =u (P i ), u 0=u (P 0),R 表示余项. 2︒ 六角网格(图14.10(b ))取P 0(x ,y )为中心,它的三个邻近节点分别为图14.9()⎪⎪⎭⎫ ⎝⎛-+-⎪⎪⎭⎫ ⎝⎛++h y h x P y h x P h y h x P 23,2,,23,2321则 R u u u h i i +∆=⎪⎭⎫⎝⎛-∑=0312334.图14.103︒ 极坐标系中的网格(图14.10(c ))取P 0(r ,θ)为中心,它的四个邻近节点分别为()()()()l r P h r P l r P h r P ++--θθθθ,,,,,,4321而拉普拉斯方程01122222=∂∂+∂∂+∂∂=θ∆u r r u r ru u的相应的差分方程为()()()011221110222134222312=⎪⎭⎫ ⎝⎛+--++++u l r h u u rh u u l r u u h 3. 抛物型方程的差分方法 考虑热传导方程的边值问题()()()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧≥==<<=><<=∂∂-∂∂0,,,,00,0,0,0,021222t t t b u t t u bx x x u t b x x u a t u μμϕ 将[0,b ]分为n 等份,每段长为∆x bn=.引两族平行线(图14.11)图14.11x =x i =i ∆x (i =0,1,2,, n )y =y j =j ∆t (j =0,1,2,, ∆t 取值见后)作成一个长方形的网格,记u (x i ,t j )为u ij ,节点(x i ,t j )为(i ,j ),在节点(i ,j )上分别用(),2,1,1,,2,1Δ2,Δ2,1,11,=-=+---++j n i x u u u t u u ji ij j i ij j i 代替22,xut u ∂∂∂∂,于是边值问题化为差分方程()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧===-===-==+----++ ,2,1,0,Δ,Δ1,,2,1,Δ,2,1,0,1,,2,10Δ2Δ21002,1,121,j t j u t j u n i x i u j n i x u u u a tu u nj j i j i ij j i ijj i μμϕ 记()22x ta ∆∆=λ,差分方程可写成 () ,2,1,1,,2,121,1,11,=-=+-+=-++j n i u u u u ji ij j i j i λλλ (1)由此可按t 增加的方向逐排求解.在第0排上u i 0的值由初值ϕ(i ∆x )确定,j +1排u i ,j +1的值可由第j 排的三点(i +1,j ),(i ,j ),(i -1,j )上的值u i +1,j , u ij ,u i -1,j 确定,而u 0,j +1,u n ,j +1已由边界条件μ1((j +1)∆t )及μ2((j +1)∆t )给定,于是可逐排计算一切节点上的u ij 值.当ϕ(x ), μ1(x )和μ2(x )充分光滑,且λ≤12时,差分方程收敛而且稳定.所以利用差分方程(1)计算时,必须使λ≤12,即()22Δ21Δx at ≤.热传导方程还可用差分方程()0Δ2Δ21,11,1,121,=+---+-++++x u u u a t u u j i j i j i ij j i 代替,此时如已知前j 排u ij 的值,为求第j +1排的u i ,j +1 必须解包含n -1个未知量u u j n j 1111,,,,+-+ 的线性代数方程组,这种差分方程称为隐式格式的差分方程,前面所提的差分方程称为显式格式差分方程.隐式格式差分方程对任意的λ都是稳定的.4. 双曲型方程的差分方法 考虑弦振动方程的第一边值问题()()()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧≥==<<=∂∂=><<=∂∂-∂∂0,,,,00),()0,(,0,0,0,02122222t t t b u t t u b x x t x u x x u t b x x u a tu μμψϕ 用矩形网格,列出对应的差分方程:()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧===-=∆=∆-==-==+--+--+-+ ,2,1,0,Δ,Δ1,,2,1),(,Δ,2,1,1,,2,1,0Δ2)(Δ22100102,1,1221,1,j t j u t j u n i x i t u u x i u j n i x u u u a t u u u nj j i i i j i ij j i j i ij j i μμψϕ 记ω=a tx∆∆与上段一样,利用u u n 022,和在第0排及第1排的已知数值(初始条件)u i 0 , u i 1可计算u i 2,然后用已知的u i 1 , u i 2及u u n 033,可计算u i 3,类似地可确定一切节点上的u ij 值.当ϕ(x ),ψ(x ),μ1(x )和μ2(x )充分光滑,且ω≤1时,差分方程收敛且稳定,所以要取∆∆t ax ≤1.二、 变分方法1. 自共轭边值问题将§3定义的共轭微分算子的概念推广到一般方程.设D 是n E 中的有界区域,S 为其边界,在D 上考虑2k 阶线性微分方程()x f x x uaLu km mi i i ni m m i i n n n=∂∂≡∑∑==++201111 ∂ 的齐次边值问题()r j u l Sj ,,2,10==式中f (x )是D 内的已知函数,l j u 是线性微分算子. 将 ⎰DvLud Ω分部积分k 次得()⎰∑⎰⎪⎪⎭⎫ ⎝⎛+=Ω=S j j j D S v R u R v u vLu d ~,Λd k 1 式中Λ(u ,v )是一个D 上的积分,其被积函数包含u ,v 的k 阶导数;R j 和 R j 是定义在边界S 上的两个线性微分算子.再将Λ(u ,v )分部积分k 次得()()⎰∑⎰⎪⎪⎭⎫⎝⎛-Ω=Λ=S k j j j D S u R v R v uL v u d ~d ,1***式中L*是一个2k 阶的微分算子,称为L 的共轭微分算子.若L=L*,则称L 为自共轭微分算子.从上面可推出格林公式()()⎰∑⎰=-=Ω-Skj jjjjDS u R v R v R u R v uL vLu 1***d ~~d 如从l j u |S =l j v |S =0可推出在边界S 上()∑==-kj jjjju R v R v R u R 1**0~~ 则称l j u |S =0为自共轭边界条件.如果微分算子及边界条件都是自共轭的,则称相应的边值问题为自共轭边值问题,此时有()0d ][=Ω-⎰DuLv vLu每个边值问题对应于某希尔伯特空间H (例如L 2(D ),见第九章§7)中的一个算子A ,其定义域M A 是H 中一线性稠密集合,它由足够次连续可微且满足边界条件的函数组成,在M A 上,Au 的数值与Lu 的数值相同,从而求解边值问题化为解算子方程Au f = 的问题.设A 为定义在实的希尔伯特空间H 中的某线性稠密集合M A 上的线性算子.若对于M A 的任意非零元素,,v u 成立(Au ,v )=(u ,Av )则称A 为对称算子.若对任意非零元素u 成立()0,>u Au 则称A 为正算子.如成立更强的不等式(Au ,u )≥r ||u ||2 (r>0)则称A 为正定算子.此处(u ,v )表示希尔伯特空间的内积,||u ||2=(u ,u ). 2. 变分原理与广义解定理 设A 是正定算子,u 是方程Au =f 在M A 上的解的充分必要条件是: u 使泛函F (u )=(Au ,u )-2(f ,u )取极小值.上述将边值问题化为等价的求泛函极值问题的方法称为能量法.在算子的定义域不够大时,泛函F (u )的极值问题可能无解.不过对于正定算子,可以开拓集合M A ,使在开拓了的集合上,泛函的极值问题有解.为开拓M A ,在M A 上引进新的内积[u ,v ]=(Au ,v ),定义模||u ||2=[u ,u ]=(Au ,u ),在模||u ||的意义下,补充极限元素,得到一个新的完备希尔伯特空间H 0,在H 0上,泛函F (u )仍然有意义,而泛函的极值问题有解.但必须注意,此时使泛函F (u )取极小的元素u 0不一定属于M A ,因此它不一定在原来的意义下满足方程Au=f 及边界条件.称u 0为广义解. 3. 极小化序列与里兹方法在处理变分问题中,极小化序列起着重要的作用.考虑泛函F (u )=(Au ,u )-2(f ,u )以d 表示泛函的极小值.设在希尔伯特空间中存在一列元素{u n } (n =1,2 ,),使()d u F n n =∞→lim则称{u n }为极小化序列.定理 若算子A 是正定的,则F (u )的每一个极小化序列既按H 空间的模也按H 0的模收敛于使泛函F (u )取极小的元素.这个定理不但指出利用极小化序列可求问题的解,而且提供一种近似解的求法,即把极小化序列中的每一个元素当作问题的近似解.设算子A 是正定的,构造极小化序列的里兹方法的主要步骤是:(1) 在线性集合M A 中选取H 0中完备的元素序列{ϕi } , (i =1,2 ,) 并要求对任意的n ,ϕ1,ϕ2,…,ϕn 线性无关.称这样的元素为坐标元素.(2) 令u a n k k k n==∑ϕ1 ,其中a k 为待定系数.代入泛函F (u ),得自变量a 1,a 2,…,a n 的函数()()()∑∑==-=nj jjn k j kjkj n f a A a a u F 11,,2,ϕϕϕ(3) 为使函数F (u n )取极小,必须()()n j a u F jn ,,2,10 ==∂∂,从而求出a k (k =1,2,…,n ).序列{u n }即为极小化序列,u n 可作为问题的近似解. 4. 里兹方法在特征值问题上的应用 算子方程Au -λu =0的非零解λ称为算子A 的特征值,对应的非零解u 称为λ所对应的特征函数. 对线性算子A ,若存在常数K ,使对任何M A 的元素ϕ成立(A ϕ,ϕ)≥K ||ϕ||2则称A 为下有界算子,正定算子是下有界的(此时K =0).记(A ϕ,ϕ)/||ϕ||2的下确界为d . 定理1 设A 为下有界对称算子,若存在不为零的元素ϕ0∈M A ,使()d A =200,ϕϕϕ则d 就是A 的最小特征值,ϕ0为对应的特征函数.于是求下有界对称算子的最小特征值问题化为变分问题,即在希尔伯特空间中求使泛函(A ϕ,ϕ)/||ϕ||2取极小的元素,或在||ϕ||=1的条件下求使泛函(A ϕ,ϕ)取极小的元素.定理2 设A 是下有界对称算子,λ1≤λ2≤…≤λn 是它的前n 个特征值,ϕ1,ϕ2,…,ϕn 是对应的标准正交特征函数,如果存在不为零的元素1+n ϕ,在附加条件(ϕ,ϕ)=1, (ϕ,ϕ1)=0, (ϕ,ϕ2)=0, …, (ϕ,ϕn )=0下使泛函(A ϕ,ϕ)取极小,则ϕn +1是算子A 的特征函数,对应的特征值()11,++=n n A ϕϕλ就是除λ1 ,,λn 外的最小的一个特征值.于是求第n +1个特征值就化为变分问题,即在附加条件(ϕ,ϕ)=1, (ϕ,ϕ1)=0, (ϕ,ϕ2)=0 ,, (ϕ,ϕn )=0 下求使泛函(A ϕ,ϕ)取极小的元素.为了利用里兹方法求特征值,在M A 中选取一列在H 0中完备的坐标元素序列{ϕi },(i =1,2 ,), 令u a n k k k n==∑ϕ1,确定a k ,使在条件 (u n ,u n )=1下,(Au n ,u n )取极小,这个问题化为求n个变元a 1,a 2,…,a n 的函数()()∑==nm k m k k m n n a a A u Au 1,,,ϕϕ在条件()()∑===nm k m k m k n n a a u u 1,1,,ϕϕ下的极值问题,一般可用拉格朗日乘数法解(见第九章§3,t ),此时()()()()()()()()()()()()0,,,,,,,,,,,,11222121111111=------n n n n n n n n n n A A A A A A ϕϕλϕϕϕϕλϕϕϕϕλϕϕϕϕλϕϕϕϕλϕϕϕϕλϕϕ的最小的根即为特征值的近似值,如果将上式的根按大小排列,就依次得后面的特征值的近似值,但精确度较差. 对一般算子方程Au -λBu=0如果A 为下有界对称算子,B 为正定算子,则()()()()()()()()()()()()0,,,,,,,,,,,,11222121111111=------n n n n n n n n n n B A B A B A B A B A B A ϕϕλϕϕϕϕλϕϕϕϕλϕϕϕϕλϕϕϕϕλϕϕϕϕλϕϕ的根就是特征值的近似值. 5. 迦辽金方法用里兹方法解数学物理问题有很多限制,最主要的限制是要求算子正定,但很多问题不一定满足这个条件,迦辽金方法弥补了这个缺陷. 迦辽金方法的主要步骤是:(1) 在M A 中选取在空间H 中完备的元素序列{ϕi } (i =1,2 ,),其中任意n 个元素线性无关,称{ϕi } (i =1,2,…)为坐标元素序列. (2) 把方程的近似解表示为u a n k k k n==∑ϕ1式中a k 是待定常数,把u n 代入方程Au=f 中的u ,两端与ϕj (j =1,2,…,n )求内积,得 a k 的n 个代数方程()()()n j f A a j n k j k k,,2,1,,1 ==∑=ϕϕϕ(3) 求出a k ,代回u n 的表达式,便得方程的近似解u n .在自共轭边值问题中,当算子是正定时,由迦辽金方法和里兹方法得到的关于a k 的代数方程组是相同的.。
数值模拟偏微分方程的三种方法介绍(有限差分方法、有限元方法、有限体积方法)I.三者简介有限差分方法(Finite Difference Methods)是数值模拟偏微分方程最早采用的方法,至今仍被广泛使用O该方法包括区域剖分和差商代替导数两个步骤。
首先将求解区域划分为差分网格,用有限个网格节点代替连续的求解区域。
其次,利用Taylor级数展开等方法将偏微分方程中的导数项在网格节点上用函数值的差商代替进行离散,从而建立以网格节点上的值为未知量的代数方程组。
该方法是一种直接将微分问题变为代数问题的近似数值解法,数学概念直观,表达简单,是发展较早且十分成熟的数值方法。
差商代替导数后的格式称为有限差分格式,从格式的精度来考虑,有一阶格式、二阶格式和高阶格式。
从差分的空间离散形式来考虑,有中心格式和迎风格式。
对于瞬态方程,考虑时间方向的离散,有显格式、隐格式、交替显隐格式等。
目前常见的差分格式,主要是以上几种格式的组合,不同的组合构成不同的差分格式。
差分方法主要适用于结构网格,网格的大小一般根据问题模型和Courant稳定条件来决定。
有限元方法(Finite Element Methods)的基础是虚位移原理和分片多项式插值。
该方法的构造过程包括以下三个步骤。
首先,利用虚位移原理得到偏微分方程的弱形式,将计算区域划分为有限个互不重叠的单元(三角形、四边形、四面体、六面体等),在每个单元上选择合适的节点作为求解函数的插值点,将偏微分方程中的变量改写成由各变量或其导数的节点值与所选用的分片插值基函数组成的线性表达式,得到微分方程的离散形式。
利用插值函数的局部支集性质及数值积分可以得到未知量的代数方程组。
有限元方法有较完善的理论基础,具有求解区域灵活(复杂区域)、单元类型灵活(适于结构网格和非结构网格)、程序代码通用(数值模拟软件多数基于有限元方法)等特点。
有限元方法最早应用于结构力学,随着计算机的发展已经渗透到计算物理、流体力学与电磁学等各个数值模拟领域。
数值计算中的偏微分方程解法偏微分方程在科学、工程和金融等领域都有广泛的应用。
在现实生活中,许多问题都涉及到偏微分方程的解法,比如天气预报、机器学习和金融衍生品定价等。
然而,解析解并不总是可行的,因此需要数值计算方法来解决这些问题。
在本文中,我们将探讨数值计算中的偏微分方程解法。
一、有限差分法有限差分法是偏微分方程数值解法中最基本的方法之一。
该方法通过将偏微分方程中的导数用差分近似公式表示出来,然后建立一个离散的空间和时间网格。
在网格上求解方程,得到数值解。
例如,考虑一个二维热传导方程:$$ \frac{\partial u}{\partial t}= \alpha \left( \frac{\partial ^2u}{\partial x^2} +\frac{\partial ^2 u}{\partial y^2} \right) $$其中,$u(x,y,t)$是温度分布,$\alpha$是热传导系数。
我们可以将该方程在空间上进行离散化,用差分近似公式表示出导数。
以二阶中心差分为例,有:$$ \frac{\partial ^2 u}{\partial x^2} \approx \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{\Delta x^2} $$$$ \frac{\partial ^2 u}{\partial y^2} \approx \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{\Delta y^2} $$其中,$u_{i,j}$表示网格点$(i,j)$处的温度。
同样地,时间上也进行离散化,用前向差分公式表示导数,即:$$ \frac{\partial u}{\partial t} \approx \frac{u_{i,j}^{n+1}-u_{i,j}^n}{\Delta t} $$将上述离散化的结果代入方程中,可以得到:$$ \frac{u_{i,j}^{n+1}-u_{i,j}^n}{\Delta t}= \alpha\left( \frac{u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n}{\Delta x^2}+\frac{u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n}{\Delta y^2} \right) $$整理得到:$$ u_{i,j}^{n+1}= u_{i,j}^n+ \frac{\alpha \Delta t}{\Delta x^2} (u_{i+1,j}^n-2u_{i,j}^n+u_{i-1,j}^n)+ \frac{\alpha \Delta t}{\Delta y^2} (u_{i,j+1}^n-2u_{i,j}^n+u_{i,j-1}^n) $$这样,我们就可以用迭代法求解上述方程,得到网格上的温度分布。
《偏微分方程数值解法》课 程 设 计题 目: 六点对称差分格式解热传导方程的初边值问题 姓 名: 王晓霜 学 院: 理学院 专 业: 信息与计算科学 班 级: 0911012 学 号:指导老师:翟方曼2012年12月14日一、题目用六点对称差分格式计算如下热传导方程的初边值问题已知其精确解为二、理论1.考虑的问题考虑一维模型热传导方程 (1.1) )(22x f xu a t u +∂∂=∂∂,T t ≤<0 其中a 为常数。
)(x f 是给定的连续函数。
(1.1)的定解问题分两类:第一,初值问题(Cauchy 问题):求足够光滑的函数()t x u ,,满足方程(1.1)和初始条件:(1.2) ()()x x u ϕ=0,, ∞<<∞-x第二,初边值问题(也称混合问题):求足够光滑的函数()t x u ,,满足方程(1.1)和初始条件:()13.1 ()()x x u ϕ=0,, l x l <<-及边值条件()23.1 ()()0,,0==t l u t u , T t ≤≤0假定()x f 和()x ϕ在相应的区域光滑,并且于()0,0,()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。
现在考虑边值问题(1.1),(1.3)的差分逼近取 N l h =为空间步长,MT =τ为时间步长,其中N ,M 是自然数, jh x x j ==, ()N j ,,1,0 =; τk y y k ==, ()M k ,,1,0 = 将矩形域G {}T t l x ≤≤≤≤=0;0分割成矩形网格。
其中 ()j i y x ,表示网格节点;h G 表示网格内点(位于开矩形G 中的网格节点)的集合;h G 表示位于闭矩形G 中的网格节点的集合;h Γ表示h G -h G 网格边界点的集合。
kj u 表示定义在网点()k i t x ,处的待求近似解,N j ≤≤0,M k ≤≤0。
注意到在节点()k i t x ,处的微商和差商之间的下列关系((,)kj k j u u x t t t ∂∂⎛⎫≡ ⎪∂∂⎝⎭): 2.区域网格剖分取空间步长lh N =和时间步长T M τ=,其中,N M 都是正整数。
用两族平行直线(0,1,,)j x x j h j N ===L 和(0,1,,)k t t k k M τ===L 将矩形域{}0;0G x l t T =≤≤≤≤分割成矩形网格,网格节点为(,)j k x t 。
以h G 表示网格内点集合,即位于开矩形的网点集合;h G 表示所有位于闭矩形G 的网点集合;h h h G G -=Γ是网格界点集合。
其次,用kj u 表示定义在网点(,)j k x t 的函数,M k N j ≤≤≤≤0,03.建立相应差分格式数值分析中,Crank-Nicolson 方法是有限差分方法中的一种,用于数值求解热方程以及形式类似的偏微分方程。
它在时间方向上是隐式的二阶方法,数值稳定。
该方法诞生于20世纪,由John Crank 与Phyllis Nicolson 发展。
①向前差分格式()24.1 ()j j j x u ϕϕ==0, k u 0=k N u =0② 向后差分格式()25.1 ()j j j x u ϕϕ==0, k u 0=k N u =0将向前差分格式和向后差分格式做算术平均,得到的差分格式称之为六点对称格式,也称为Grank-Nicholson 格式:()26.1 ()j j j x u ϕϕ==0, k u 0=k N u =0进一步,()36.1 2r-11++kj u +()r +11+k j u 2r -11+-k j u =2r k j u 1++()r -1k j u 2r +k j u 1-+j f τ 按层计算:首先,取0=k ,则利用初值()j j j x u ϕϕ==0和边值k u 0=k N u =0,来确定出第一层的1j u ,1,,1,0-=N j ,即求解方程组:2r -11+j u +()r +11j u 2r -11-j u =2r 01+j u +()r -10j u 2r +01-j u +j f τ 1,,1,0-=N j ,k u 0=k N u =0。
求出1j u ,在由()18.1',取1=k ,可利用1j u ,解出2j u ,1,,1,0-=N j 。
如此下去,即可逐层算出所有kj u ,1,,1,0-=M k 。
若记三、截断误差()=u R kj ()()[]k j k j h Lu t x u L -,3=()22h O +τ。
注意:21+⎪⎪⎭⎫ ⎝⎛∂∂k j t u =()()τk j k j t x u t x u ,,1-+()2τO + 又 两式相加而21+⎪⎪⎭⎫ ⎝⎛∂∂k j t u 2122+⎪⎪⎭⎫ ⎝⎛∂∂⋅=k j x u a +j f 故有()22h O ++τ。
四、稳定性与收敛性抛物方程的两层差分格式可以统一写成向量形式:(1.7) 1k k AU BU F τ+=+其中1111(,), (,,)k k k T N N U u u F f f --==,A 和B 是1N -阶矩阵。
我们假定A 可逆,即(1.7)是唯一可解的。
对于显格式,A 等于单位矩阵I 。
三层格式可以通过引入新变量1k kk U W U -⎛⎫= ⎪⎝⎭化成两层格式。
假设差分解的初始值(其实可以是任一层的值)有误差,以后各层计算没有误差,让我们来考察初始误差对以后各层的影响。
令{}k U 和{}k V 分别是以0U 和0V 为初始值由差分格式(1.7)得到的两组差分解,则k k k W U V =-满足(1.8) 1k k AW BW += 因此,按初值稳定应该意味着10k W K W +≤。
这就导致如下定义:假设0F =,我们称差分格式(2.1)按初值稳定,如果存在正常数0τ和K ,使得以下不等式成立:(1.8) 10k U K U +≤, 010, 0, 0N U R k T τττ-∀∈≤≤<≤ 这里∙是1N R -上的某一个范数,例如类似地,假设00U =,我们称差分格式(2.1)按右端稳定,如果存在正常数0τ和K ,使得以下不等式成立:(1.8) 1k U K F +≤, 00, 0k T τττ∀≤≤<≤可以证明,差分格式若按初值稳定,则一定按右端稳定。
因此,这时我们简单地称差分格式稳定。
前面讨论的向前差分格式①当网比12r ≤时稳定,当12r >时不稳定。
这就意味着给定空间步长h 以后,时间步长τ必须足够小,才能保证稳定。
而向后差分格式②和Grank-Nicholson 格式(1.6)则对任何网比r 都是稳定的,时间步长τ可以取得大一些,从而提高运算效率。
如果某个差分格式的截断误差当τ和h 趋于0时随之趋于0,则称这个差分格式是相容的。
可以证明:若差分格式是相容的和稳定的,则它是收敛的,并且差分解与微分解之间误差的阶等于截断误差的阶。
因此,对任何网比,向后差分格式(1.6)有收敛阶2()O h τ+。
五、结论对于扩散方程(包括许多其他方程),可以证明Crank-Nicolson 方法无条件稳定。
但是,如果时间步长与空间步长平方的比值过大(一般地,大于1/2),近似解中将存在虚假的振荡或衰减。
基于这个原因,当要求大时间步或高空间分辨率的时候,往往会采用数值精确较差的后向欧拉方法进行计算,这样即可以保证稳定,又避免了解的伪振荡。
由本题可以总结出,抛物型方程的六点对称差分法所得的数值解能够较好地逼近方程的精确 解,且区域剖分得越细,即步长越小,数值解与精确解的误差就越小,数值解越逼近精确解。
六、附录 取51=h ,1001=τ,则212a r h τ==,满足稳定性条件 取110h =,1400τ=,则212a r h τ==,满足稳定性条件 另取120h =,11600τ=,则212a r h τ==,亦满足稳定性条件 程序代码format longa=2;l=1;T=1;N=5;M=100;h=l/N;to=T/M;r=(a*to)/h^2;for j=1:N+1x(j)=(j-1)*h;for k=1:M+1t(k)=(k-1)*to;u(j,k)=exp(x(j)+2*t(k));endend u %求解精确解for j=1:N+1x(j)=(j-1)*h;us(j,1)=exp(x(j));endfor k=1:M+1t(k)=(k-1)*to;us(1,k)=exp(2*t(k));us(N+1,k)=exp(1+2*t(k));endfor k=2:M+1for j=2:Nus(j,k)=r*us(j-1,k-1)+(1-2*r)*u(j,k-1)+r*us(j+1,k-1);endendus %求解数值解for k=1:M+1for j=1:N+1R(j,k)=abs(u(j,k)-us(j,k));endendR %计算误差Rmax=max(max(R)) %求误差的最大值精确解与数值解的比较:x=0:0.1:1;hold onplot(x,u(:,M+1),'b');plot(x,us(:,M+1),'y');title('t=1,h=1/10,τ=1/400时精确解和数值解的比较') text(0.05,21,'蓝:精确解');text(0.05,20,'黄:数值解');hold off取不同步长时的误差比较:x=0:1/5:1;y=0:1/10:1;z=0:1/20:1;hold onplot(x,R(:,M+1),'b');hold offN分别取5, 10,20①当N=5时,误差最大值:②当N=10时,误差最大值:③当N=20时,误差最大值:。