ADI(交替方向隐格式)求解二维抛物方程(含matlab程序)

合集下载

二维双曲型方程的隐格式解法

二维双曲型方程的隐格式解法

本科毕业论文(设计)题目二维双曲线型方程的交替方向隐格式解法院(系)数学系专业数学及应用数学学生姓名周玲玲学号 10022156指导教师陈淼超职称讲师论文字数 9500完成日期: 2014年6月8日巢湖学院本科毕业论文(设计)诚信承诺书本人郑重声明:所呈交的本科毕业论文(设计),是本人在导师的指导下,独立进行研究工作所取得的成果。

除文中已经注明引用的内容外,本论文不含任何其他个人或集体已经发表或撰写过的作品成果。

对本文的研究做出重要贡献的个人和集体,均已在文中以明确方式标明。

本人完全意识到本声明的法律结果由本人承担。

本人签名:日期:巢湖学院本科毕业论文 (设计)使用授权说明本人完全了解巢湖学院有关收集、保留和使用毕业论文 (设计)的规定,即:本科生在校期间进行毕业论文(设计)工作的知识产权单位属巢湖学院。

学校根据需要,有权保留并向国家有关部门或机构送交论文的复印件和电子版,允许毕业论文 (设计)被查阅和借阅;学校可以将毕业论文(设计)的全部或部分内容编入有关数据库进行检索,可以采用影印、缩印或扫描等复制手段保存、汇编毕业,并且本人电子文档和纸质论文的内容相一致。

保密的毕业论文(设计)在解密后遵守此规定。

本人签名:日期:导师签名:日期:二维双曲型方程的交替方向隐格式解法摘要在解决偏微分方程中二维双曲线型方程的问题求解初边值问题时,显示格式的稳定性是有条件的,并且多维的稳定性条件更严,为得到稳定性好的格式,隐式格式受到重视。

用隐式格式求解二维问题得到的线性方程组其系数矩阵为宽带状,因此求解不甚便利,采用交替方向隐式(ADI )格式可以避免此问题。

本文以基本概念和基本方法为主,同时结合算例实现算法。

第一部分介绍二维双曲线型方程的交替方向隐格式解法的基本概念,引入本文的研究对象——二维波动方程: ,x R ∈,],0(T t ∈第二部分介绍上述方程的二维双曲线型方程的交替方向隐格式及这种格式的存在性、收敛性及稳定性。

用向后euler法求三维抛物型方程的adi差分格式

用向后euler法求三维抛物型方程的adi差分格式

文章标题:向后Euler法在三维抛物型方程中的应用及ADI差分格式探讨在数学和科学计算领域,求解三维抛物型方程是一个重要且复杂的问题。

本文将从数值方法的角度出发,讨论如何使用向后Euler法来求解三维抛物型方程,并探讨ADI差分格式在该过程中的应用。

1. 三维抛物型方程的数值求解在数学建模和科学计算中,三维抛物型方程的数值求解是一个广泛应用的问题。

三维抛物型方程一般具有形如$\frac{\partial u}{\partial t} = \nabla \cdot (D\nabla u) + f$的形式,其中$D$为扩散系数,$f$为源项函数。

由于方程复杂性,传统的解析方法难以得到精确解,因此需要借助数值方法来进行求解。

2. 向后Euler法向后Euler法是一种常用的数值方法,用于离散化时间导数。

其基本思想是将时间导数用差分近似替代,通过迭代计算得到时间步数上的解。

对于三维抛物型方程,可以将时间方向的偏导数用向后Euler法进行离散化,从而得到数值解。

3. ADI差分格式ADI(Alternating Direction Implicit)差分格式是一种常用的隐式差分方法,用于求解多维偏微分方程。

其核心思想是将多维偏微分方程拆分为一维方程的求解,通过交替方向隐式差分得到整体方程的数值解。

在三维抛物型方程的求解中,ADI差分格式能够有效地提高计算效率和数值稳定性。

4. 主题回顾与总结通过本文的介绍,我们了解了向后Euler法在三维抛物型方程求解中的重要性和应用。

ADI差分格式作为一种高效的数值方法,为复杂方程的求解提供了可行的途径。

对于三维抛物型方程,我们可以利用向后Euler法结合ADI差分格式,得到高质量、深度和广度兼具的数值解。

5. 个人观点与理解在数值计算中,选择合适的数值方法对于求解复杂方程至关重要。

向后Euler法作为一种简单而有效的数值方法,为我们提供了一种直观且可行的思路。

而ADI差分格式则在多维问题的求解中发挥着重要作用,其交替方向求解的思想能够有效地提高计算效率和数值稳定性。

六点对称法,ADI法,预校法,和LOD法解二维抛物线方程

六点对称法,ADI法,预校法,和LOD法解二维抛物线方程

GAGGAGAGGAFFFFAFAF偏微分數值解法實驗報告實驗名稱:六點對稱格式,ADI 法,預校法,LOD 法解二維拋物線方程的初值問題實驗成員: 吳興 楊敏 姚榮華 于瀟龍 余凡 鄭永亮實驗日期:2013年5月17日 指導老師:張建松一、 實驗內容用六點對稱格式,ADI 法,預校法和LOD 法求解二維拋物線方程的初值問題:21(),(,)(0,1)(0,1),0,4(0,,)(1,,)0,01,0,(,0,)(,1,)0,01,0(,,0)sin cos .xx yy y y u u u x y G t t u y t u y t y t u x t u x t x t u x y x y ππ∂⎧=+∈=⨯>⎪∂⎪⎪==<<>⎨⎪==<<>⎪=⎪⎩已知(精確解為:2(,,)sin cos exp()8u x y t x y t πππ=-)設(0,1,,),(0,1,,),(0,1,,)j k n x jh j J y kh k K t n n N τ======差分解為GAGGAGAGGAFFFFAFAF,nj k u ,則邊值條件為:0,,,0,1,1,0,0,1,,,,0,1,,n n k J k nn n nj j j K j K u u k K u u u u j J-⎧===⎪⎨===⎪⎩初值條件為:0,sin cos j k j k u x y ππ=取空間步長12140h h h ===,時間步長11600τ=網比。

1: ADI 法:由第n 层到第n+1层计算分成两步:先先第n 層到n+1/2層,對uxx 用向后差分逼近,對uyy 用向前差分逼近,對uyy 用向后差分逼近,于是得到了如下格式:11112222,,1,,1,,1,,1221222,,2-22=21()n n n n n n n nj kj kj kj k j kj k j k j k n nx j k y j k hh hτδδ+++++-+-+-+-+=+uu uuuu u u (+) (1)u u1111111222,,1,,1,,1,,12212212,,2-22=21()n n n n n n n n j kj kj kj k j kj k j k j k n n x j k y j k hh hτδδ++++++++-+-++-+-+=+u u uuuu u u (+) (2)u u其中j,k=1,2,…,M-1,n=0,1,2,…,上标n+1/2表示在t=tn+1/2+(n+1/2)τ取值。

二维抛物型方程的交替方向隐格式

二维抛物型方程的交替方向隐格式

二维抛物型方程的交替方向隐格式二维抛物型方程的交替方向隐格式是一种数值解法,用于求解二维抛物型方程的数值解。

这种解法将二维问题分解为两个一维问题,并采用隐式差分方法来求解。

具体来说,二维抛物型方程可以表示为:u_t = a^2 u_{xx} + b^2 u_{yy} + f(x,y,u)其中,u是待求解的函数,t是时间,a和b是两个不同的参数,f是右侧的非线性函数。

为了求解这个问题,我们可以采用交替方向隐式差分方法,将问题分解为两个一维问题:1、在x方向上,从左到右扫描每一行数据,更新每个点的u值。

这个过程可以使用隐式差分方法来实现:u^[i,j]^(n+1) = u^[i,j]^(n) + dt a^2 (u^[i+1,j]^(n) - 2u^[i,j]^(n) + u^[i-1,j]^(n)) / h^2 + dt f^[i,j]^(n) 其中,u^[i,j]^(n)表示第n个时间步中,位置为(x[i],y[j])的点的u值,h是x方向上的步长,dt是时间步长,f^[i,j]^(n)表示第n个时间步中,位置为(x[i],y[j])的点上的非线性函数f的值。

2. 在y方向上,从上到下扫描每一列数据,更新每个点的u值。

这个过程也可以使用隐式差分方法来实现:u^[i,j]^(n+1) = u^[i,j]^(n) + dt b^2 (u^[i,j+1]^(n) - 2u^[i,j]^(n) + u^[i,j-1]^(n)) / h^2 + dt f^[i,j]^(n) 其中,u^[i,j]^(n)表示第n个时间步中,位置为(x[i],y[j])的点的u值,h是y方向上的步长,dt是时间步长,f^[i,j]^(n)表示第n个时间步中,位置为(x[i],y[j])的点上的非线性函数f的值。

通过交替更新每个点的u值,我们可以逐步逼近方程的数值解。

这种解法具有二阶精度和稳定性,可以应用于多种二维抛物型方程的问题。

二维抛物方程的有限差分法

二维抛物方程的有限差分法

二维抛物方程的有限差分法二维抛物方程的有限差分法摘要二维抛物方程是一类有广泛应用的偏微分方程,由于大部分抛物方程都难以求得解析解,故考虑采用数值方法求解。

有限差分法是最简单又极为重要的解微分方程的数值方法。

本文介绍了二维抛物方程的有限差分法。

首先,简单介绍了抛物方程的应用背景,解抛物方程的常见数值方法,有限差分法的产生背景和发展应用。

讨论了抛物方程的有限差分法建立的基础,并介绍了有限差分方法的收敛性和稳定性。

其次,介绍了几种常用的差分格式,有古典显式格式、古典隐式格式、Crank-Nicolson隐式格式、Douglas差分格式、加权六点隐式格式、交替方向隐式格式等,重点介绍了古典显式格式和交替方向隐式格式。

进行了格式的推导,分析了格式的收敛性、稳定性。

并以热传导方程为数值算例,运用差分方法求解。

通过数值算例,得出古典显式格式计算起来较简单,但稳定性条件较苛刻;而交替方向隐式格式无条件稳定。

关键词:二维抛物方程;有限差分法;古典显式格式;交替方向隐式格式FINITE DIFFERENCE METHOD FORTWO-DIMENSIONAL PARABOLICEQUATIONAbstractTwo-dimensional parabolic equation is a widely used class of partial differential equations. Because this kind of equation is so complex, we consider numerical methods instead of obtaining analytical solutions. finite difference method is the most simple and extremely important numerical methods for differential equations. The paper introduces the finite difference method for two-dimensional parabolic equation.Firstly, this paper introduces the background and common numerical methods for Parabolic Equation, Background and development of applications. Discusses the basement for the establishment of the finite difference method for parabolic equation And describes the convergence and stability for finite difference method.Secondly, Introduces some of the more common simple differential format,for example, the classical explicit scheme, the classical implicit scheme, Crank-Nicolson implicit scheme, Douglas difference scheme, weighted six implicit scheme and the alternating direction implicit format. The paper focuses on the classical explicit scheme and the alternating direction implicit format. The paper takes discusses the derivation convergence,and stability of the format . The paper takes And the heat conduction equation for the numerical example, using the differential method to solve. Through numerical examples, the classical explicit scheme is relatively simple for calculation, with more stringent stability conditions; and alternating direction implicit scheme is unconditionally stable.Keywords:Two-dimensional Parabolic Equation; Finite-Difference Method; Eclassical Explicit Scheme; Alternating Direction Implicit Scheme1绪论1.1课题背景抛物方程是一类特殊的偏微分方程,二维抛物方程的一般形式为u Lu t ∂=∂ (1-1)其中1212((,,))((,,))(,,)(,,)(,,)u u u u u u L a x y t a x y t b x y t b x y t C x y t x x y y x y∂∂∂∂∂∂=++++∂∂∂∂∂∂ 120,0,0a a C >>≥。

ADI(交替方向隐格式)求解二维抛物方程(含matlab程序)

ADI(交替方向隐格式)求解二维抛物方程(含matlab程序)

ADI 法求解二维抛物方程学校:中国石油大学(华东) 学院:理学院 姓名:张道德 时间:2013.4.271、ADI 法介绍作为模型,考虑二维热传导方程的边值问题:(3.6.1),0,,0(,,0)(,)(0,,)(,,)(,0,)(,,)0t xx yy u u u x y l t u x y x y u y t u l y t u x t u x l t ϕ=+<<>⎧⎪=⎨⎪====⎩取空间步长1hM,时间步长0,作两族平行于坐标轴的网线:,,,0,1,,,j k x x jh y y kh j k M =====将区域0,x y l ≤≤分割成2M 个小矩形。

第一个ADI 算法(交替方向隐格式)是Peaceman 和Rachford (1955)提出的。

方法:由第n 层到第n+1层计算分为两步:(1) 第一步: 12,12n j k xx yy u +从n->n+,求u 对向后差分,u 向前差分,构造出差分格式为:1(3.6.1)11112222,,1,,1,,1,,1221222,,2-22=21()n n n n n n n n j kj kj kj k j kj k j k j k n n x j k y j k hhhτδδ+++++-+-+-+-+=+uu uuuu u u (+)u u(2) 第二步:12,12n j k xx yy u +从n+->n+1,求u 对向前差分,u 向后差分,构造出差分格式为:2(3.6.1)1111111222,,1,,1,,1,,12212212,,2-22=21()n n n n n n n n j kj kj kj k j kj k j k j k n n x j k y j k hh hτδδ++++++++-+-++-+-+=+uu uuuu u u (+)u u其中1211,1,,1,0,1,2,,()22n j k M n n n τ+=-=+=+上表表示在t=t 取值。

二维抛物型方程的交替方向隐格式

二维抛物型方程的交替方向隐格式

二维抛物型方程的交替方向隐格式在数学领域中,二维抛物型方程是一类重要的偏微分方程,它们在众多实际问题的数学建模中起着关键作用。

对于这类方程,交替方向隐格式是一种常用且有效的数值求解方法。

本文将详细介绍二维抛物型方程及其交替方向隐格式的原理和应用,希望能为读者提供一份生动、全面且有指导意义的参考材料。

首先,我们来了解什么是二维抛物型方程。

通常,二维抛物型方程可以表示为以下形式:∂u/∂t = α(∂²u/∂x² + ∂²u/∂y²) + βu + f(x, y, t)其中,u是未知函数,t是时间变量,x和y是空间变量,α和β是常数,f(x, y, t)是已知函数。

二维抛物型方程广泛应用于物理、工程、生物等领域的问题求解,比如热传导、扩散、扩散反应等。

为了求解这类方程,数学家们开发了各种求解方法,交替方向隐格式是其中一种。

交替方向隐格式是一种时间和空间交替迭代的求解方法,它通过将二维抛物型方程离散化为一组代数方程,然后通过迭代求解这些代数方程得到数值解。

具体来说,交替方向隐格式先将时间方向离散化,将时间变量t划分为一系列离散时间步长。

然后,对于每个时间步长,交替方向隐格式将二维抛物型方程中的时间导数∂u/∂t进行近似。

最常用的近似方法是向后差分格式,即用u(n+1) - u(n)来近似∂u/∂t,其中u(n)表示第n个时间步长的数值解,u(n+1)表示第n+1个时间步长的数值解。

这样,二维抛物型方程可以离散为一组代数方程。

接下来,交替方向隐格式将空间方向离散化,将空间变量x和y划分为一系列离散网格点。

然后,在空间离散化的基础上,通过引入交替方向(例如,先按x方向更新,再按y方向更新,或者反之)和隐格式(例如,使用向后差分格式近似二阶导数项),将二维抛物型方程中的空间导数进行近似。

通过交替迭代求解这组离散代数方程,我们可以得到二维抛物型方程的数值解。

当离散网格点的数量足够多时,数值解将趋近于方程的解。

matlab求解二维抛物线型偏微分方程

matlab求解二维抛物线型偏微分方程

一、简介二维抛物线型偏微分方程是一类常见的偏微分方程,在科学与工程领域有着重要的应用。

利用Matlab求解二维抛物线型偏微分方程是一种常见的数值求解方法,它可以帮助我们快速地得到方程的数值解,并对问题进行分析和研究。

二、二维抛物线型偏微分方程的一般形式二维抛物线型偏微分方程一般可表示为:∂u/∂t = ∂^2u/∂x^2 + ∂^2u/∂y^2 + f(x, y, t)其中,u是未知函数,f(x, y, t)是给定的函数,代表外力或源项。

这类偏微分方程描述了许多现实世界中的问题,如传热传质、扩散反应等。

三、使用Matlab求解二维抛物线型偏微分方程的基本步骤1. 网格划分:将求解区域进行离散化,构建网格。

2. 离散化方程:将偏微分方程进行差分处理,得到一个离散的代数方程组。

3. 求解代数方程组:利用Matlab中的求解器求解得到问题的数值解。

4. 后处理:对数值解进行可视化和分析,得出结论并进行讨论。

四、具体例子考虑二维热传导方程:∂u/∂t = α(∂^2u/∂x^2 + ∂^2u/∂y^2)其中,α是热传导系数。

假设我们要求解一个长方形区域上的热传导问题,边界条件已知,初值条件也已给出。

我们可以利用Matlab 进行数值求解。

五、Matlab代码示例下面是一个简单的Matlab代码示例,用于求解二维热传导方程:定义问题的基本参数Lx = 1; 区域长度Ly = 1; 区域宽度Nx = 100; 网格数Ny = 100;dx = Lx / Nx; 网格步长dy = Ly / Ny;alpha = 0.01; 热传导系数dt = 0.001; 时间步长Nt = 1000; 时间步数初始化温度场u = zeros(Nx, Ny);设置边界条件和初值条件...用有限差分方法离散化方程for n = 1:Nt计算下一个时间步的温度场...end可视化结果...后处理...六、结论利用Matlab求解二维抛物线型偏微分方程是一种高效、便捷的数值方法,能够帮助我们快速地得到问题的数值解,并对问题进行分析和研究。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

ADI 法求解二维抛物方程学校:中国石油大学(华东) 学院:理学院 姓名:张道德 时间:2013.4.271、ADI 法介绍作为模型,考虑二维热传导方程的边值问题:(3.6.1),0,,0(,,0)(,)(0,,)(,,)(,0,)(,,)0t xx yy u u u x y l t u x y x y u y t u l y t u x t u x l t ϕ=+<<>⎧⎪=⎨⎪====⎩取空间步长1hM,时间步长0,作两族平行于坐标轴的网线:,,,0,1,,,j k x x jh y y kh j k M =====将区域0,x y l ≤≤分割成2M 个小矩形。

第一个ADI 算法(交替方向隐格式)是Peaceman 和Rachford (1955)提出的。

方法:由第n 层到第n+1层计算分为两步:(1) 第一步: 12,12n j k xx yy u +从n->n+,求u 对向后差分,u 向前差分,构造出差分格式为:1(3.6.1)11112222,,1,,1,,1,,1221222,,2-22=21()n n n n n n n n j kj kj kj k j kj k j k j k n n x j k y j k hhhτδδ+++++-+-+-+-+=+uu uuuu u u (+)u u(2) 第二步:12,12n j k xx yy u +从n+->n+1,求u 对向前差分,u 向后差分,构造出差分格式为:2(3.6.1)1111111222,,1,,1,,1,,12212212,,2-22=21()n n n n n n n n j kj kj kj k j kj k j k j k n n x j k y j k hh hτδδ++++++++-+-++-+-+=+uu uuuu u u (+)u u其中1211,1,,1,0,1,2,,()22n j k M n n n τ+=-=+=+上表表示在t=t 取值。

假定第n 层的,n j ku已求得,则由1(3.6.1)求出12,n j k u +,这只需按行(1,,1)j M =-解一些具有三对角系数矩阵的方程组;再由2(3.6.1)求出1,n j k u +,这只需按列(1,,1)k M =-解一些具有三对角系数矩阵的方程组,所以计算时容易实现的。

2、数值例子(1)问题用ADI 法求解二维抛物方程的初边值问题:21(),(,)(0,1)(0,1),0,4(0,,)(1,,)0,01,0,(,0,)(,1,)0,01,0(,,0)sin cos .xx yy y y u u u x y G t t u y t u y t y t u x t u x t x t u x y x y ππ∂⎧=+∈=⨯>⎪∂⎪⎪==<<>⎨⎪==<<>⎪=⎪⎩ 已知(精确解为:2(,,)sin cos exp()8u x y t x y t πππ=-)设(0,1,,),(0,1,,),(0,1,,)j k n x jh j J y kh k K t n n N τ======差分解为,n j k u ,则边值条件为:0,,,0,1,1,0,0,1,,,,0,1,,n nk J k nn n nj j j K j K u u k Ku u u u j J-⎧===⎪⎨===⎪⎩初值条件为:0,sin cos j k j k u x y ππ=取空间步长12140h h h ===,时间步长11600τ=网比21r h τ==。

用ADI 法分别计算到时间层1t =。

(2)计算过程从n 到n+1时,根据边值条件:0,,0,0,1,,n n k J k u u k K ===,已经知道第0列和第K 列数值全为0。

(1)12,12n j k xx yy u +从n->n+,求u 对向后差分,u 向前差分,构造出差分格式为:11112222,,1,,1,,1,,1221222,,2-221=1621()16n n n n n n n n j kj kj kj k j kj k j k j k n n x j k y j k hhhτδδ+++++-+-+-+-+=+u u u u uu u u (+)u u从而得到:1112221,,1,,1,,1111111(1)(1)321632321632n n n n n n j k j k j k j k j k j k r r r r r r ++++-+--++-=+--u u u u u u ,其中1,2,,1,1,2,,1j J k K =-=-即按行用追赶法求解一系列下面的三对角方程组:121,122,123,123,122,121,(1)(1)111163211113216321111321632111132163211113216321113216n k n k n k n J k n J k n J k J J r r r r r r rr r r r r r r r r ++++-+-+--⨯-⎡⎤⎡+-⎢⎥⎢⎢⎥⎢⎢⎥-+-⎢⎢⎥⎢⎢⎥⎢⎢⎥-+-⎢⎢⎥⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥⎢⎥-+⎣⎢⎥⎣⎦u u u u u u 123321(1)1(1)1J J J J J f f f f f f ----⨯-⨯⎤⎥⎥⎡⎤⎥⎢⎥⎥⎢⎥⎥⎢⎥⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎢⎥⎦ 又根据边值条件得:,0,1,1,,,0,1,,nnnnj j j K j K u u u u j J -===,解出第0行,0nj u 和第K 行,,(0,1,,)n j K u j J =。

(2)第二步:12,12n j k xx yy u +从n+->n+1,求u 对向前差分,u 向后差分,构造出差分格式为:1111111222,,1,,1,,1,,12212212,,2-221=1621()16n n n n n n n n j kj kj kj k j kj k j k j k n n x j k y j k hhhτδδ++++++++-+-++-+-+=+uu u uuu u u (+)u u从而得到:111111222,1,,11,,1,111111(1)(1)321632321632n n n n n n j k j k j k j k j k j k r r r r r r +++++++-+--++-=+--u u u u u u ,其中1,2,,1,1,2,,1j J k K =-=-又根据边值条件得:,0,1,1,,,0,1,,nnnnj j j K j K u u u u j J -===,从而得到:,0,1,1,0n nj j n nj K j K u u u u -⎧-=⎪⎨-+=⎪⎩其中(0,1,,)j J = 即按列用追赶法求解一系列下面的三对角方程组:1,01,11,21,31,31,21111113216321111321632111132163211113216321111321632111132163211n j n j n j n j n j K n j K K Kr r r r r r r r r r r rr r r r r r +++++-+-⨯-⎡⎤⎢⎥⎢⎥-+-⎢⎥⎢⎥-+-⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥-+-⎢⎥⎢⎥⎢⎥-+-⎢⎥⎢⎥-⎢⎥⎣⎦u u u u u u u 12343211,111,1K K n K j K K K n j K K f f f f f f f f --+--⨯+⨯⎡⎤⎢⎥⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦u (3) 求解结果从而得到误差的范数为:1- 范数:0.233770443573713; 2-范数:0.196807761631447; ∞-范数:0.327253314506086(3.4)图像(3.4.1)数值解图像:(3.4.2)精确解图像:(5)主要程序(5.1)主程序%*************************************************************%main_chapter主函数%信息10-2——张道德%学号:10071223clccleara = 0; b=1; %x取值范围c=0; d=1; %y取值范围tfinal = 1; %最终时刻t=1/1600;%时间步长;h=1/40;%空间步长r=t/h^2;%网比x=a:h:b;y=c:h:d;%**************************************************************%精确解m=40;u1=zeros(m+1,m+1);for i=1:m+1,for j=1:m+1u1(j,i) = uexact(x(i),y(j),1);endend%数值解u=ADI(a,b,c,d,t,h,tfinal);%***************************************************************** %绘制图像figure(1) ;mesh(x,y,u1)figure(2); mesh(x,y,u)%误差分析error=u-u1;norm1=norm(error,1);norm2=norm(error,2);norm00=norm(error,inf);%*****************************************************************(5.2)ADI函数%****************************************************************% 用ADI法求解二维抛物方程的初边值问题% u_t = 1/16(u_{xx} + u_{yy})(0,1)*(0,1)% 精确解: u(t,x,y) = sin(pi*x) sin(pi*y)exp(-pi*pi*t/8) %****************************************************************function [u]=ADI(a,b,c,d,t,h,tfinal )%(a , b) x取值范围%(c, d) y取值范围%tfinal最终时刻%t时间步长;%h空间步长r=t/h^2;%网比m=(b-a)/h;%n=tfinal/t; %x=a:h:b;y=c:h:d;%****************************************************************** %初始条件u=zeros(m+1,m+1);for i=1:m+1,for j=1:m+1u(j,i) = uexact(x(i),y(j),0);endend%******************************************************************u2=zeros(m+1,m+1);a=-1/32*r*ones(1,m-2);b=(1+r/16)*ones(1,m-1);aa=-1/32*r*ones(1,m);cc=aa;aa(m)=-1;cc(1)=-1;bb=(1+r/16)*ones(1,m+1);bb(1)=1;bb(m+1)=1;for i=1:n%************************************************************** %从n->n+1/2,u_{xx}向后差分,u_{yy}向前差分for j=2:mfor k=2:md(k-1)=1/32*r*(u(j,k+1)-2*u(j,k)+u(j,k-1))+u(j,k);end% 修正第一项与最后一项,但由于第一项与最后一项均为零,可以省略%d(1)=d(1)+u1(j,1);d(m-1)=d(m-1)+u1(j,m+1);u2(j,2:m)=zhuiganfa(a,b,a,d);endu2(1,:)=u2(2,:);u2(m+1,:)=u2(m,:);%************************************************************** %从n->n+1,u_{xx}向前差分,u_{yy}向后差分for k=2:mdd(1)=0;dd(m+1)=0;for j=2:mdd(j)=1/32*r*(u2(j+1,k)-2*u2(j,k)+u2(j-1,k))+u2(j,k);endu(:,k)=zhuiganfa(aa,bb,cc,dd);end%**************************************************************** u2=u;end%********************************************************************(5.3)“追赶法”程序%******************************************************************** %追赶法function [x]=zhuiganfa(a,b,c,d)%对角线下方的元素,个数比A少一个% %对角线元素%对角线上方的元素,个数比A少一个%d为方程常数项%用追赶法解三对角矩阵方程r=size(a);m=r(2);r=size(b);n=r(2);if size(a)~=size(c)|m~=n-1|size(b)~=size(d)error('变量不匹配,检查变量输入情况!');end%%%LU分解u(1)=b(1);for i=2:nl(i-1)=a(i-1)/u(i-1);u(i)=b(i)-l(i-1)*c(i-1);v(i-1)=(b(i)-u(i))/l(i-1);end%追赶法实现%%%求解Ly=d,"追"的过程y(1)=d(1);for i=2:ny(i)=d(i)-l(i-1)*y(i-1);end%%%求解Ux=y,"赶"的过程x(n)=y(n)/u(n);for i=n-1:-1:1x(i)=y(i)/u(i);x(i)=(y(i)-c(i)*x(i+1))/u(i);end%********************************************************************(5.4)精确解函数%t时刻,u的取值;function [ f]=uexact(x,y,t)f=sin(x*pi)*cos(y*pi)*exp(-pi*pi/8*t);%********************************************************************。

相关文档
最新文档