ADI(交替方向隐格式)求解二维抛物方程(含matlab程序)
二维双曲型方程的隐格式解法
本科毕业论文（设计）题目二维双曲线型方程的交替方向隐格式解法院（系）数学系专业数学及应用数学学生姓名周玲玲学号 10022156指导教师陈淼超职称讲师论文字数 9500完成日期: 2014年6月8日巢湖学院本科毕业论文(设计)诚信承诺书本人郑重声明：所呈交的本科毕业论文(设计)，是本人在导师的指导下，独立进行研究工作所取得的成果。
除文中已经注明引用的内容外，本论文不含任何其他个人或集体已经发表或撰写过的作品成果。
对本文的研究做出重要贡献的个人和集体，均已在文中以明确方式标明。
本人完全意识到本声明的法律结果由本人承担。
本人签名：日期：巢湖学院本科毕业论文 (设计)使用授权说明本人完全了解巢湖学院有关收集、保留和使用毕业论文 (设计)的规定，即：本科生在校期间进行毕业论文(设计)工作的知识产权单位属巢湖学院。
学校根据需要，有权保留并向国家有关部门或机构送交论文的复印件和电子版，允许毕业论文 (设计)被查阅和借阅；学校可以将毕业论文(设计)的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存、汇编毕业，并且本人电子文档和纸质论文的内容相一致。
保密的毕业论文(设计)在解密后遵守此规定。
本人签名：日期：导师签名：日期：二维双曲型方程的交替方向隐格式解法摘要在解决偏微分方程中二维双曲线型方程的问题求解初边值问题时，显示格式的稳定性是有条件的，并且多维的稳定性条件更严，为得到稳定性好的格式，隐式格式受到重视。
用隐式格式求解二维问题得到的线性方程组其系数矩阵为宽带状，因此求解不甚便利，采用交替方向隐式（ADI ）格式可以避免此问题。
本文以基本概念和基本方法为主，同时结合算例实现算法。
第一部分介绍二维双曲线型方程的交替方向隐格式解法的基本概念，引入本文的研究对象——二维波动方程： ，x R ∈，],0(T t ∈第二部分介绍上述方程的二维双曲线型方程的交替方向隐格式及这种格式的存在性、收敛性及稳定性。
用向后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法解二维抛物线方程
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 法求解二维抛物方程学校：中国石油大学（华东） 学院：理学院 姓名：张道德 时间：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求解二维抛物线型偏微分方程是一种常见的数值求解方法，它可以帮助我们快速地得到方程的数值解，并对问题进行分析和研究。
二、二维抛物线型偏微分方程的一般形式二维抛物线型偏微分方程一般可表示为：∂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求解二维抛物线型偏微分方程是一种高效、便捷的数值方法，能够帮助我们快速地得到问题的数值解，并对问题进行分析和研究。
二维抛物线方程数值解法(ADI隐式交替法)方法
ADI隐式交替法三种解法及误差分析（一般的教材上只说第一种）理论部分参看孙志忠：偏微分方程数值解法注意：1.最好不要直接看程序，中间很多公式很烦人的（一定要小心），我写了两天，终于写对了。
2.中间：例如r*(u(i-1,m1,k)+u(i+1,m1,k))形式写成分形式：r*u(i-1,m1,k)+r*u(i+1,m1,k)后面会出错，我也不是很清楚为什么，可能由于舍入误差，或者大数吃掉小数的影响。
3.下面有三个程序4.具体理论看书，先仔细看书（孙志忠：偏微分方程数值解法）或者网上搜一些理论。
Matlab程序：1.function [u u0 p e x y t]=ADI1(h1,h2,m1,m2,n)%ADI解二维抛物线型偏微分方程(P-R交替隐式，截断)%此程序用的是追赶法解线性方程组%h1为空间步长，h2为时间步长%m1,m2分别为x方向，y方向网格数，n为时间网格数%p为精确解，u为数值解,e为误差%定义u0(i,j,k)=u(i,j,k+1/2),因为矩阵中，i,j,k必须全为整数x=(0:m1)*h1+0;%定义x0,y0,t0是为了f(x,t)~=0的情况%y=(0:m2)*h1+0;t=(0:n)*h2+0; t0=(0:n)*h2+1/2*h2;for k=1:n+1for i=1:m2+1for j=1:m1+1f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k));endendendfor i=1:m2+1for j=1:m1+1u(i,j,1)=exp(0.5*(x(j)+y(i)));endendfor k=1:n+1for i=1:m2+1u(i,[1 m1+1],k)=[exp(0.5*y(i)-t(k)) exp(0.5*(1+y(i))-t(k))]; u0(i,[1 m1+1],k)=[exp(0.5*y(i)-t0(k)) exp(0.5*(1+y(i))-t0(k))] ;endendfor k=1:n+1for j=1:m1+1u([1 m2+1],j,k)=[exp(0.5*x(j)-t(k)) exp(0.5*(1+x(j))-t(k))]; u0([1 m2+1],j,k)=[exp(0.5*x(j)-t0(k)) exp(0.5*(1+x(j))-t0(k))];endendr=h2/(h1*h1);r1=2*(1-r);r2=2*(1+r);for k=1:n %外循环，先固定每一时间层，每一时间层上解一线性方程组% for i=2:m2a=-r*ones(1,m1-1);c=a;a(1)=0;c(m1-1)=0;b=r2*ones(1,m1-1);d(1)=r*u0(i,1,k)+r*(u(i-1,2,k)+u(i+1,2,k))+r1*u(i,2,k)+...h2*f(i,2,k);for l=2:m1-2d(l)=r*(u(i-1,l+1,k)+u(i+1,l+1,k))+r1*u(i,l+1,k)+...h2*f(i,l+1,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m1-1)=r*u0(i,m1+1,k)+r*(u(i-1,m1,k)+u(i+1,m1,k))...+r1*u(i,m1,k)+h2*f(i,m1,k);for l=1:m1-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu0(i,m1,k)=d(m1-1)/b(m1-1); %回代过程%for l=m1-2:-1:1u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k))/b(l);endendfor j=2:m1a=-r*ones(1,m2-1);c=a;a(1)=0;c(m2-1)=0;b=r2*ones(1,m2-1);d(1)=r*u(1,j,k+1)+r*(u0(2,j-1,k)+u0(2,j+1,k))+r1*u0(2,j,k)+...h2*f(2,j,k);for l=2:m2-2d(l)=r*(u0(l+1,j-1,k)+u0(l+1,j+1,k))+r1*u0(l+1,j,k)+... h2*f(l+1,j,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m2-1)=r*u(m2+1,j,k+1)+r*(u0(m2,j-1,k)+u0(m2,j+1,k))...+r1*u0(m2,j,k)+h2*f(m2,j,k);for l=1:m2-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu(m2,j,k+1)=d(m2-1)/b(m2-1); %回代过程%for l=m2-2:-1:1u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);endendendfor k=1:n+1for i=1:m2+1for j=1:m1+1p(i,j,k)=exp(0.5*(x(j)+y(i))-t(k)); %p为精确解e(i,j,k)=abs(u(i,j,k)-p(i,j,k)); %e为误差endendend2.function [u p e x y t]=ADI2(h1,h2,m1,m2,n)%ADI解二维抛物线型偏微分方程(D'Yakonov交替方向隐格式)%此程序用的是追赶法解线性方程组%h1为空间步长，h2为时间步长%m1,m2分别为x方向，y方向网格数，n为时间网格数%p为精确解，u为数值解,e为误差%定义u0(i,j,k)=u'(i,j,k)(引入的过渡层),因为矩阵中，i,j,k必须全为整数x=(0:m1)*h1+0;y=(0:m2)*h1+0;t=(0:n)*h2+0;t0=(0:n)*h2+1/2*h2;%定义t0是为了f(x,y,t)~=0的情况% for k=1:n+1for i=1:m2+1for j=1:m1+1f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k));%编程时-t0(k)写成了+t0(k),导致错误;endendend%初始条件for i=1:m2+1for j=1:m1+1u(i,j,1)=exp(0.5*(x(j)+y(i)));endend%边界条件for k=1:n+1for i=1:m2+1u(i,[1 m1+1],k)=[exp(0.5*y(i)-t(k)) exp(0.5*(1+y(i))-t(k))];endendr=h2/(h1*h1);r4=1+r;r5=r/2;for k=1:nfor i=2:m2u0(i,[1 m1+1],k)=r4*u(i,[1 m1+1],k+1)-r5*(u(i-1,[1 m1+1],... k+1)+u(i+1,[1 m1+1],k+1));endendfor k=1:n+1for j=1:m1+1u([1 m2+1],j,k)=[exp(0.5*x(j)-t(k)) exp(0.5*(1+x(j))-t(k))];endendr1=r-r*r;r2=2*(r-1)*(r-1);r3=r*r/2;for k=1:n %外循环，先固定每一时间层，每一时间层上解一线性方程组% for i=2:m2a=-r*ones(1,m1-1);c=a;a(1)=0;c(m1-1)=0;b=2*r4*ones(1,m1-1);d(1)=r*u0(i,1,k)+r1*(u(i-1,2,k)+u(i,1,k)+u(i+1,2,k)+...u(i,3,k))+r2*u(i,2,k)+r3*(u(i-1,1,k)+...u(i+1,1,k)+u(i-1,3,k)+u(i+1,3,k))+2*h2*f(i,2,k);for l=2:m1-2d(l)=r1*(u(i-1,l+1,k)+u(i,l,k)+u(i+1,l+1,k)+...u(i,l+2,k))+r2*u(i,l+1,k)+r3*(u(i-1,l,k)+...u(i+1,l,k)+u(i-1,l+2,k)+u(i+1,l+2,k))+2*h2*f(i,l+1,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m1-1)=r*u0(i,m1+1,k)+r1*(u(i-1,m1,k)+u(i,m1-1,k)+...u(i+1,m1,k)+u(i,m1+1,k))+r2*u(i,m1,k)+...r3*(u(i-1,m1-1,k)+...u(i+1,m1-1,k)+u(i-1,m1+1,k)+u(i+1,m1+1,k))+2*h2*f(i,m1,k);for l=1:m1-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);end%回代过程%u0(i,m1,k)=d(m1-1)/b(m1-1);for l=m1-2:-1:1u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k))/b(l);endendfor j=2:m1a=-r*ones(1,m2-1);c=a;a(1)=0;c(m2-1)=0;b=2*r4*ones(1,m2-1);d(1)=r*u(1,j,k+1)+2*u0(2,j,k);for l=2:m2-2d(l)=2*u0(l+1,j,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m2-1)=2*u0(m2,j,k)+r*u(m2+1,j,k+1);for l=1:m2-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu(m2,j,k+1)=d(m2-1)/b(m2-1); %回代过程%for l=m2-2:-1:1u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);endendendfor k=1:n+1for i=1:m2+1for j=1:m1+1p(i,j,k)=exp(0.5*(x(j)+y(i))-t(k)); %p为精确解e(i,j,k)=abs(u(i,j,k)-p(i,j,k)); %e为误差endendend3.function [u u0 p e x y t]=ADI5(h1,h2,m1,m2,n)%ADI解二维抛物线型偏微分方程(P-R交替隐式，未截断)%此程序用的是追赶法解线性方程组%h1为空间步长，h2为时间步长%m1,m2分别为x方向，y方向网格数，n为时间网格数%p为精确解，u为数值解,e为误差%定义u0(i,j,k)=u(i,j,k+1/2),因为矩阵中，i,j,k必须全为整数x=(0:m1)*h1+0;%定义x0,y0,t0是为了f(x,t)~=0的情况%y=(0:m2)*h1+0;t=(0:n)*h2+0; t0=(0:n)*h2+1/2*h2;for k=1:n+1for i=1:m2+1for j=1:m1+1f(i,j,k)=-1.5*exp(0.5*(x(j)+y(i))-t0(k));endendendfor i=1:m2+1for j=1:m1+1u(i,j,1)=exp(0.5*(x(j)+y(i)));endendfor k=1:n+1for i=1:m2+1u(i,[1 m1+1],k)=[exp(0.5*y(i)-t(k)) exp(0.5*(1+y(i))-t(k))]; u1(i,[1 m1+1],k)=[exp(0.5*y(i)-t0(k)) exp(0.5*(1+y(i))-t0(k))] ;endendr=h2/(h1*h1);r1=2*(1-r);r2=r/4;r3=2*(1+r);for k=1:nfor i=2:m2u0(i,[1 m1+1],k)=u1(i,[1 m1+1],k)-r2*(u(i-1,[1 m1+1],k+1)-...2*u(i,[1 m1+1],k+1)+u(i+1,[1 m1+1],k+1)-u(i-1,[1 m1+1],k)+...2*u(i,[1 m1+1],k)-u(i+1,[1 m1+1],k));endendfor k=1:n+1for j=1:m1+1u([1 m2+1],j,k)=[exp(0.5*x(j)-t(k)) exp(0.5*(1+x(j))-t(k))];endendfor k=1:n %外循环，先固定每一时间层，每一时间层上解一线性方程组%for i=2:m2a=-r*ones(1,m1-1);c=a;a(1)=0;c(m1-1)=0;b=r3*ones(1,m1-1);d(1)=r*u0(i,1,k)+r*(u(i-1,2,k)+u(i+1,2,k))+r1*u(i,2,k)+... h2*f(i,2,k);for l=2:m1-2d(l)=r*(u(i-1,l+1,k)+u(i+1,l+1,k))+r1*u(i,l+1,k)+...h2*f(i,l+1,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m1-1)=r*u0(i,m1+1,k)+r*(u(i-1,m1,k)+u(i+1,m1,k))...+r1*u(i,m1,k)+h2*f(i,m1,k);for l=1:m1-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu0(i,m1,k)=d(m1-1)/b(m1-1); %回代过程%for l=m1-2:-1:1u0(i,l+1,k)=(d(l)-c(l)*u0(i,l+2,k))/b(l);endendfor j=2:m1a=-r*ones(1,m2-1);c=a;a(1)=0;c(m2-1)=0;b=r3*ones(1,m2-1);d(1)=r*u(1,j,k+1)+r*(u0(2,j-1,k)+u0(2,j+1,k))+r1*u0(2,j,k)+...h2*f(2,j,k);for l=2:m2-2d(l)=r*(u0(l+1,j-1,k)+u0(l+1,j+1,k))+r1*u0(l+1,j,k)+... h2*f(l+1,j,k);%输入部分系数矩阵，为0的矩阵元素不输入%一定要注意输入元素的正确性endd(m2-1)=r*u(m2+1,j,k+1)+r*(u0(m2,j-1,k)+u0(m2,j+1,k))...+r1*u0(m2,j,k)+h2*f(m2,j,k);for l=1:m2-2 %开始解线性方程组消元过程a(l+1)=-a(l+1)/b(l);b(l+1)=b(l+1)+a(l+1)*c(l);d(l+1)=d(l+1)+a(l+1)*d(l);endu(m2,j,k+1)=d(m2-1)/b(m2-1); %回代过程%for l=m2-2:-1:1u(l+1,j,k+1)=(d(l)-c(l)*u(l+2,j,k+1))/b(l);endendendfor k=1:n+1for i=1:m2+1for j=1:m1+1p(i,j,k)=exp(0.5*(x(j)+y(i))-t(k)); %p为精确解 e(i,j,k)=abs(u(i,j,k)-p(i,j,k)); %e为误差endendend[up e x y t]=ADI2(0.01,0.001,100,100,1000);surf(x,y,e(:,:,1001)) t=1的误差曲面下面是三种方法的误差比较：1.[u u0 p e x y t]=ADI1(0.1,0.1,10,10,10)((P-R交替隐式，截断)截断中间过渡层用u(i,j,k+1/2）代替）（t=1时的误差）2.[u u0 p e x y t]=ADI5(0.1,0.1,10,10,10)(P-R交替隐式，未截断)(未截断过渡层u(i,j,)’=u(i,j,k+1/2)-h2^2/4*dy^2dtu(i,j,k+1/2);)3.[u p e x y t]=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隐格式) surf(x,y,e(:,:,11))(表示t=1时的误差）下面是相关数据：1: [u u0 p e x y t]=ADI1(0.1,0.1,10,10,10)((P-R交替隐式，截断)截断中间过渡层用u(i,j,k+1/2)代替）e(:,:,11) =Columns 1 through 60 0 0 0 0 00 0.00040947 0.00025182 0.00019077 0.00017112 0.000176040 0.00057359 0.00042971 0.00035402 0.00032565 0.000336280 0.00066236 0.00054689 0.00047408 0.00044596 0.000462670 0.00072152 0.00062001 0.00055081 0.00052442 0.000545530 0.00076164 0.0006576 0.00058522 0.00055732 0.000579840 0.00078336 0.00065993 0.00057557 0.00054161 0.000562090 0.00078161 0.00061872 0.00051646 0.00047429 0.000489640 0.00073621 0.0005148 0.00039979 0.00035439 0.000363130 0.00056964 0.00031688 0.00022051 0.0001884 0.000191920 0 0 0 0 02.[u u0 p e x y t]=ADI5(0.1,0.1,10,10,10)(P-R交替隐式，未截断)(未截断过渡层u(i,j,)’=u(i,j,k+1/2)-h2^2/4*dy^2dtu(i,j,k+1/2);)e(:,:,11) =Columns 1 through 60 0 0 0 0 00 0.00027006 0.00016305 0.00012104 0.0001071 0.000109950 0.00037754 0.00027817 0.0002253 0.00020483 0.000211160 0.00043539 0.00035386 0.00030207 0.00028124 0.00029140 0.00047398 0.00040104 0.00035113 0.00033111 0.000344050 0.0005003 0.00042535 0.00037309 0.0003519 0.000365710 0.00051479 0.00042699 0.00036681 0.00034164 0.00035410 0.00051415 0.00040056 0.00032887 0.0002985 0.000307640 0.00048504 0.0003335 0.00025411 0.0002221 0.000227060 0.00037609 0.00020532 0.00013956 0.00011718 0.000119020 0 0 0 0 03.[u p e x y t]=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隐格式)e(:,:,11) =Columns 1 through 60 0 0 0 0 00 8.6469e-006 1.4412e-005 1.8364e-005 2.091e-005 2.2174e-0050 1.4412e-005 2.4777e-005 3.2047e-005 3.6716e-005 3.8961e-0050 1.8364e-005 3.2047e-005 4.1789e-005 4.8054e-005 5.1008e-0050 2.091e-005 3.6716e-005 4.8054e-005 5.5353e-005 5.8764e-0050 2.2174e-005 3.8961e-005 5.1008e-005 5.8764e-005 6.2389e-0050 2.2118e-005 3.8698e-005 5.0523e-005 5.8126e-005 6.171e-0050 2.055e-005 3.5581e-005 4.6157e-005 5.2942e-005 5.6197e-0050 1.707e-005 2.8951e-005 3.7128e-005 4.2365e-005 4.4952e-0050 1.0851e-005 1.7698e-005 2.2265e-005 2.5203e-005 2.672e-0050 0 0 0 0 01.[u u0 p e x y t]=ADI1(0.1,0.1,10,10,10)((P-R交替隐式，截断)截断中间过渡层用u(i,j,k+1/2)代替）Columns 7 through 110 0 0 0 00.00020348 0.00026228 0.00038338 0.00066008 00.00038607 0.00048321 0.00064717 0.00091668 00.00052635 0.00064203 0.00081637 0.0010517 00.0006174 0.00074272 0.00092111 0.0011417 00.00065651 0.00078964 0.00097724 0.0012051 00.00064051 0.00078116 0.00098594 0.0012433 00.00056474 0.00070822 0.00093332 0.0012478 00.00042547 0.00055526 0.00078616 0.0011844 00.00022735 0.00030946 0.00049004 0.00092402 00 0 0 0 02.[u u0 p e x y t]=ADI5(0.1,0.1,10,10,10)(P-R交替隐式，未截断)(未截断过渡层u(i,j,)’=u(i,j,k+1/2)-h2^2/4*dy^2dtu(i,j,k+1/2);)Columns 7 through 110 0 0 0 00.00012826 0.00016798 0.00024986 0.00043637 00.00024444 0.00031023 0.00042173 0.00060513 00.00033401 0.00041257 0.00053179 0.00069358 00.00039216 0.00047742 0.00059986 0.00075263 00.00041704 0.00050761 0.00063642 0.00079439 00.00040657 0.0005021 0.00064226 0.00081984 00.00035784 0.000455 0.00060828 0.00082334 00.00026866 0.00035628 0.00051263 0.0007824 00.00014262 0.00019789 0.00031956 0.0006113 00 0 0 0 0 3.[u p e x y t]=ADI2(0.1,0.1,10,10,10)(D'Yakonov交替方向隐格式) Columns 7 through 110 0 0 0 02.2118e-005 2.055e-005 1.707e-005 1.0851e-005 03.8698e-005 3.5581e-005 2.8951e-005 1.7698e-005 05.0523e-005 4.6157e-005 3.7128e-005 2.2265e-005 05.8126e-005 5.2942e-005 4.2365e-005 2.5203e-005 06.171e-005 5.6197e-005 4.4952e-005 2.672e-005 06.1116e-005 5.5803e-005 4.4817e-005 2.6785e-005 05.5803e-005 5.1239e-005 4.1529e-005 2.5153e-005 04.4817e-005 4.1529e-005 3.42e-005 2.126e-005 02.6785e-005 2.5153e-005 2.126e-005 1.3869e-005 00 0 0 0 0。
二维抛物线方程数值解法方法
二维抛物线方程数值解法方法二维抛物线方程是数学中常见的偏微分方程之一，表示一个二维区域中的传热、扩散、波动等现象。
求解二维抛物线方程通常需要使用数值解法，其中ADI（Alternating Direction Implicit）隐式交替法是一种有效的数值解法之一、本文将详细介绍ADI隐式交替法的原理和步骤。
首先，我们来看一下二维抛物线方程的一般形式：∂u/∂t=α(∂²u/∂x²+∂²u/∂y²)+f(x,y,t)其中u(x,y,t)是未知函数，表示二维区域中的一些物理量（如温度），α是扩散系数，f(x,y,t)是源（即外力项）。
ADI隐式交替法的基本思想是将偏微分方程中的时间导数项进行分离，然后通过交替方向的半隐式差分逼近来求解。
具体步骤如下：1.将时间导数项进行分离，得到两个互相独立的方程：∂u/∂t=α(∂²u/∂x²)+α(∂²u/∂y²)+f(x,y,t)等价于：∂u/∂t=α(∂²u/∂x²)+f(x,y,t)，考虑到y方向上的导数已经被分离出来。
∂u/∂t=α(∂²u/∂y²)+f(x,y,t)，考虑到x方向上的导数已经被分离出来。
2. 对第一个方程应用隐式差分（比如Crank-Nicolson差分格式），得到一个关于未知函数u的线性方程组Ax=b。
3.利用ADI思想，对第二个方程同样进行隐式差分，得到另一个关于未知函数u的线性方程组Ay=c。
4.对线性方程组Ax=b和Ay=c进行求解，得到u的一个时间层的近似解。
5.重复步骤2至4，直到得到整个时间区间内的数值解。
ADI隐式交替法的关键在于递归地求解两个互相独立的线性方程组。
在每个时间步长中，先求解一个方向（比如x方向）上的方程组，然后用该方程组的解代入另一个方向（比如y方向）的方程组中求解。
通过交替进行，可以逼近二维抛物线方程的解。
李荣华二维抛物方程的初边值问题matlab
李荣华二维抛物方程的初边值问题在数学和工程领域中扮演着重要角色。
通过使用Matlab，我们能够更深入地理解和解决这一问题。
在本文中，我们将从基础概念开始，逐步深入讨论李荣华二维抛物方程的初边值问题，并结合Matlab进行实际分析和解决。
一、初边值问题的概念和涵义1. 初边值问题的定义初边值问题是指在偏微分方程中，除了部分边界条件外，还需给出一些初始条件。
李荣华二维抛物方程的初边值问题即是在方程中给定初值条件和边界条件的问题。
2. 李荣华二维抛物方程简介李荣华二维抛物方程是描述热传导、扩散等现象的数学模型，在物理学和工程领域有着广泛的应用。
它的形式通常为一个关于未知函数和其在空间和时间上的导数的方程。
3. Matlab在初边值问题中的应用Matlab作为一种强大的数学建模和仿真工具，能够帮助我们求解各种类型的偏微分方程，包括初边值问题。
通过Matlab，我们可以对李荣华二维抛物方程进行数值求解和模拟，从而更清晰地理解问题本质。
二、李荣华二维抛物方程的初边值问题解析1. 方程的建立我们需要建立李荣华二维抛物方程的数学模型，并给定相应的初值条件和边界条件。
在Matlab中，我们可以通过定义矩阵和方程的形式来实现这一过程。
2. 数值求解我们利用Matlab提供的数值求解方法，如有限差分法或有限元法，对初边值问题进行求解。
通过程序的编写和运行，我们可以得到方程的数值解，并对物理过程有更直观的认识。
3. 结果分析在求得数值解之后，我们可以对结果进行分析和可视化。
通过Matlab 的绘图功能，我们能够直观地观察李荣华二维抛物方程的演化过程，以及初值条件和边界条件对解的影响。
三、个人观点和总结在本文中，我们深入探讨了李荣华二维抛物方程的初边值问题，并结合Matlab进行了实际分析和解决。
通过对问题的全面评估和数值求解，我们对方程的特性和行为有了更深入的理解。
个人观点上，我认为Matlab是一个非常强大的工具，能够帮助我们更直观、高效地处理数学问题。
基于matlab解抛物型方程的交替隐方向p-r差分格式的实现
基于matlab解抛物型方程的交替隐方向p-r差分格式的实现1. 引言1.1 概述本文旨在利用MATLAB中的抛物型方程解析方法，具体实现交替隐方向p-r差分格式。
抛物型方程是一类常见的偏微分方程，在科学计算和工程领域中有着广泛的应用。
该类方程描述了许多自然界和社会系统中的动态过程，如热传导、扩散、弹性形变等。
而交替隐方向p-r差分格式则是一种高效解法，适用于求解抛物型方程。
1.2 文章结构本文将按以下结构展开详细论述：- 第2节将简要介绍抛物型方程及其解析方法概述，并特别关注MATLAB在此过程中的应用。
- 第3节将深入探讨交替隐方向p-r差分格式的原理，并对其稳定性和精确度进行分析。
- 第4节将重点阐述基于MATLAB实现交替隐方向p-r差分格式的步骤，包括空间离散化方法选择与实现、时间离散化方法选择与实现、以及迭代求解过程描述与收敛性分析。
- 最后，第5节将呈现数值实验设置，并展示数值结果，同时对结果进行讨论。
1.3 目的本文的目的在于通过MATLAB解析抛物型方程，并实现交替隐方向p-r差分格式，从而提供一种高效、稳定、精确的数值计算方法。
此研究对于处理抛物型方程相关问题具有实际应用意义，为科学计算和工程领域中的相关研究提供了指导和借鉴。
我们期望该研究能够拓展数值计算方面的知识，促进在实践中解决复杂系统动态过程模拟与分析的能力。
2. 抛物型方程解析2.1 抛物型方程简介抛物型方程是一类常见的偏微分方程，它描述了许多自然现象和数学模型中的动态行为。
一般而言，抛物型方程包括一个时间变量和多个空间变量，并且通常具有二阶时间导数和二阶或更高阶的空间导数。
典型的抛物型方程包括热传导方程、扩散方程和波动方程等。
2.2 解析方法概述解析方法是指通过使用数学分析和解析推导来求解偏微分方程的方法。
在抛物型方程的解析研究中，常用的方法包括变量分离法、相似变量法、格林函数法等。
这些方法基于物理建模和数学推导，可以得到精确的解或者近似解。
ADI(交替隐式求求解抛物型偏微分方程)程序
% 右端函数: f(t,x,y) = 1/2*exp( 1/2*(x+y) - t ) %
% 初始条件: g1 = exp( 1/2*(x+y) ) %
h = (b-a)/n; dt=h;
h1 = h*h;
r=dt/h1;
x=a:h:b; y=c:h:d;
%-- Initial condition:
t = 0;
for i=1:m+1,
for j=1:m+1,
u1(i,j) = uexact(t,x(i),y(j));
% %
% Files needed for the test: %
% %
% %
% Example of ADI Method for 2D heat equation %
% %
% u_tt = u_{xx} + u_{yy} + f(x,t) %
% a<x<b;c<y<d %
% Test problme: %
end %-- Finished with the loop in time
%----------- Data analysis ----------------------------------
for i=1:m+1,
for j=1:n+1,
ue(i,j) = uexact(tfinal,x(i),y(j));
e = max(max(abs(u1-ue))) % The infinity error.
一类二维抛物型方程的ADI格式
一类二维抛物型方程的ADI格式【摘要】本文针对一类二维抛物型方程，建立了一个在空间和时间方向上均具有二阶精度的ADI格式，并分析其稳定性. 比较以往算法，此格式具有精度较高，无条件稳定等优点，同时，该方法通过求解两个线性代数方程得到原问题的解，避免了非线性迭代运算，提高了计算效率.【关键词】二维抛物型方程；ADI格式；稳定性；截断误差1.引言抛物型偏微分方程在研究热传导过程、一些扩散现象及电磁场传播等许多问题中都有广泛的应用，对这一类方程数值解法的研究一直是科研工作者关注的热点问题之一，其中高精度的有限差分方法更是受到了越来越多的重视. 考虑如下的初边值问题[1]：其中，为常数.在文献[1]中对问题（1）建立了差分格式，格式的截断误差阶为.本文将在文献[1]的基础上进一步研究问题（1）的高效差分格式，建立了一个高精度的交替方向隐式差分格式（即ADI格式），提高了时间方向上的精度，并给出相应的稳定性分析。
2.差分格式的建立为了得到问题（1）的有限差分格式，首先将求解区域进行网格剖分，结点为. 选取正整数L和N，并令为方向上的网格步长，为方向上的网格步长，记假定第层的已知，则由第（Ⅰ）步，通过解一个三对角线性代数方程组求出，再由第（Ⅱ）步，再解一个三对角线性代数方程组即可求出. 所以，只要利用追赶法求解两个三对角线性代数方程组即可，此时计算量与储存量都较小. 另外，在处理边界条件时，为了提高精度，采取中心差分，这样会出现虚拟值，此时，只要再与格式中的方程联立，即可消去虚拟值[2].3. 稳定性分析下面采用von Newmann方法[3]对上述D格式进行稳定性分析. 一般地，低阶项不影响差分格式的稳定性，所以不妨略去项，并对（3）、（5）式消去中间变量得：利用Taylor展开式求误差，可知此处建立的D格式的截断误差阶为.参考文献：[1]管秋琴.一类二维抛物型方程的有限差分格式[J]. 赤峰学院学报（自然科学版）. 2010，26（1）：7.[2]Richtmyer R D ，Morton KW. Difference methods for initial - value problems （2nd edition）[J ] . John wiley &sons. 1967.[3]戴嘉尊，邱建贤. 微分方程数值解法[M]. 南京：东南大学出版社.2002.；安庆师范学院青年科研基金项目（KJ201020）。
二维变系数抛物型方程的一个高阶ADI差分格式
二维变系数抛物型方程的一个高阶ADI差分格式马小霞;颜晓琳;陈汝栋【摘要】A high accuracy alternation direction implicit scheme (ADI) for solving the two-dimensional parabolic equations is presented, and the scheme is absolutely stable and the truncation error is O(τ2+h4). The experiments show the scheme is effective and advantage, and the theory is right by a numerical example.%针对二维变系数抛物型方程，构造出了一个高精度、恒稳定的交替方向隐式（ADI）差分格式，格式的截断误差阶达O（τ2+h4）。
通过数值实验，验证了理论分析的正确性和差分格式的精确性与有效性。
【期刊名称】《天津工业大学学报》【年(卷),期】2014(000)001【总页数】5页(P77-80,84)【关键词】抛物型方程;ADI格式;截断误差;恒稳定【作者】马小霞;颜晓琳;陈汝栋【作者单位】焦作大学基础部，河南焦作 454003;天津工业大学理学院，天津300387;天津工业大学理学院，天津 300387【正文语种】中文【中图分类】O241.82抛物型方程在处理废料污染、渗透、驱动、海水入侵以及半导体等工程实际问题中有着广泛的应用，因此研究其高精度、高稳定和计算量较小的数值解法具有重要的意义.用有限差分方法研究这类问题的数值方法目前已做了许多工作[1-5].但这些工作大多是对常系数而言的.文献[4]中对二维变系数抛物型方程数值方法仅对系数依赖于一个变量的情况进行了研究，本文的研究是对系数依赖于两个变量的情形进行的.应用Taylor展开、算子方法[6]以及粘结系数法[7]得到了一个高精度（截断误差阶达O（τ2+h4））、恒稳定的ADI格式.格式的建立和稳定性分析都比文献[4]简单得多，文末的数值实验证明了本文理论分析的正确性和所得格式的精确性与有效性.考虑如下的二维变系数非齐次抛物型方程初边值问题其中：0＜c1≤a（x，y）≤c2；Γ为Ω的边界.设τ=Δt=T/N为时间步长，h=Δx=Δy=1/M为空间步长，N、M均为正整数.表示在节点（jh， kh，nτ）处的网函数值，微分方程问题（1）—（3）的解函数为u（x，y，t），并记u（jh，kh，nτ）=u（j，k，n），由Taylor展开式在（j，k，n+1）处考虑方程（1），根据粘结系数法[7]有其中：分别为关于x、y的二阶中心差分算子；η1、η2、η3均为待定参数.适当选取这些参数使得式（9）—式（11）成立.由于即式（9）—式（11）的3个式子均成立，由式（8）、（14）、（16）和（17）得：将所得的等式左端加上右端加上，经过简单的计算可知，所加项的误差阶均不超过O（τ2+h4），则可得在（18）中略去小量项，并令，则可建立如下差分格式：易知差分格式（19）—（21）的截断误差为对格式（19）—（21）引进中间过渡量，可以得到问题（1）—（3）的ADI格式：当第n层上的值已知时，对于固定的k，（22）是关于的三对角线性方程组，可用追赶法求解.在解方程组（22）时，需要边值和（0≤k≤M-1），这可从方程（23）得到，在方程（23）中，分别令j=0和j=M，有当过渡层已求得时，对于固定的j，式（23）是关于的三对角线性方程组，也可用追赶法求解.解方程组（23）所需要的边界值和（0≤k≤M-1）可直接从已知的边界值条件得到，即：为对格式（22）—（23）进行稳定性分析，可将过渡层变量消去，即得格式（19），故两格式有相同的截断误差和稳定性质.因为差分格式关于初值稳定一定可以推出关于右端稳定，故只对齐次格式进行讨论，根据稳定性分析的Fourier方法[21]，令将上式代入格式（19）式的齐次形式中，经计算整理，并利用关系式这里，可得格式（19）的传播因子为由s1的取值范围可知，所以|G（s1，s2）|根据Von neumann条件知格式（19）也即格式（22）—（23）恒稳定.综合上节与本节的论述，并根据Lax的稳定性与收敛性等价原理可得：定理本文的ADI格式恒稳定且以o（τ2+h4）的收敛阶收敛.在区域D：{0≤x，y≤1，t≥0}上对初值边问题用本文格式求数值解并与精确解u（x，y，t）=e-2tsin（x+y），相比较，取h=1/10，τ=rh2=r/100，r=1/2计算到n=200时的结果如表1所示.由表1结果看出，对所取的r，本文格式解与精确解均有很好的吻合，这与理论分析完全一致.本文构造了二维变系数抛物型方程的ADI格式，数值例子表明，它具有易于计算（能三对角追赶法求解）、精度高（较现有格式）、无条件稳定等特点.本文构造出的格式可以推广到任意维变系数抛物型方程，对当前关于变系数抛物型方程差分格式的构造与研究有重要的理论和实践意义.【相关文献】[1]DAI Weizhong，NASSAR Raja.A Second-order ADI scheme for three-dimensional parabolic differential equations[J].Numer Math，1964（6）：428-453.[2]DAI Weizhong，NASSAR pact ADI method for solving parabolic differential equations[J].Numer Methods Partial Differential Eq，2002（18）：129-142.[3] 孙志忠，李雪玲.反应扩散方程的紧交替方向差分格式[J].计算数学，2005，19（2）：209-224.[4] 孙志忠，李雪玲.二维变系数反应扩散方程的紧交替方向差分格式[J].高等学校计算数学学报，2006，28（1）：83-95.[5] 马菊意，杨辉.二维抛物型方程的一个高精度PC格式[J].工程数学学报，2008，25（2）：373-376.[6]余德浩，汤华中.微分方程数值解法[M].北京：科学出版社，2003：102-174.[7] 戴嘉尊，邱建贤.微分方程数值解法[M].南京：东南大学出版社，2003：91-93[8]李立康，於崇华，朱政华.微分方程数值方法[M].上海：复旦大学出版社，1999：208-213.。
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 取值。
matlab报告——用matlab研究抛体运动
用matlab研究抛体运动2. 用matlab研究抛体运动引论MATLAB语言是一种集数值计算、符号运算、可视化建模、仿真和图形处置等多种功能的高级语言。
利用MATLAB模拟物理现象为咱们解决问题提供了一种新的方式，利用其方便的数值计算和作图功能，能够方便的模拟一些物理进程。
关于处置非线性问题，既能进行数值求解，又能绘制有关曲线，方便有效，基于其功能壮大，界面友善，语言自然，交互性强等优势，已成为教学和科研中最基础的软件之一，利用其解决复杂的数值计算问题，能够减少工作量，节约时刻，图形绘制问题，真实直观，能够加深明白得，提高工作效率将物体以必然的初速度向空中抛出，仅在重力作用下物体所作的运动，它的初速度不为零，可分为平抛运动和斜抛运动。
物理上提出的“抛体运动”是一种理想化的模型，即把物体看成质点，抛出后只考虑重力作用，忽略空气阻力。
抛体运动加速度恒为重力加速度，相等的时刻内速度转变量相等，而且速度转变的方向始终是竖直向下的。
抛体运动及应用、实验设计思路一、理论分析一样的处置方式是将其分解为水平方向和竖直方向，平抛运动水平方向是匀速直线运动，竖直方向是自由落体运动，斜抛运动水平方向是匀速直线运动，竖直方向是竖直上抛运动，在任意方向上分解有正交分解和非正交分解两种情加速度及位移等进行相应分析。
不管如何分解，都必需把运动的独立性和独立作用原理结合进行系统分解，即将初速度、受力情、加速度及位移等进行相应分析。
斜抛运动： 水平方向速度αcos 0v v x= （1）竖直方向速度gt v v y -=αsin 0 （2）水平方向位移 tx v αcos 0= （3）竖直方向位移 2021cos gtt y v -=α （4）平抛运动： 水平方向速度v v x 0=（5）竖直方向速度gt v y = (6) 水平方向位移tx v 0= (7)竖直方向位移221gtv y = (8)合速度t g vv vv y xt 42202241+=+=（9）合速度方向与水平夹角β：v v v gt tg x y 0==β (10)合位移yxs 22+=(11)位移方向与水平夹角α：02v gttg ss xy==α (12)设某一抛射体的初速度为0v ，抛射角为θ，将其运动在X,Y 轴上进行正交分解，水平方向速度0cos x v v θ=(13)竖直方向0sin y v v gt θ=- (14) 质点的坐标(,)x y 是0()cos()x t t v θ= (15)201()sin 2y t t gt v θ=- (16)从上两式消去t ，便得质点的轨迹运动方程2220tan 2cos gx y x v θθ=-t (17)抛射体能达到的最大高度为220sin 2H g vθ= (18)其抵达最大高度所需时刻为0sin T gv θ= (19)空中飞行时刻为0sin 22t T g v θ== (20)抛射体的最大射程为2sin 2X gv θ= (21)它跟初速度0v 和抛射角θ有关，在抛射角θ不变的情形下，射程x 与20v 成正比，因此射程随初速度的增大而增大。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 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);%********************************************************************。