抛物型方程隐格式的程序实现
一种大时间步长求解抛物型方程的显—隐式格式的数值试验
一种大时间步长求解抛物型方程的显—隐式格式的数值试验向新华;史汉生
【期刊名称】《教学与研究（南京）》
【年(卷),期】2003(024)001
【摘要】并行计算中，在用显-隐式格式求解抛物型方程时，由于受显式格式的控制，网比比较小，导致时间步长也比较小，不能充分显示出隐式格式的优势，本文提出一种新的算法，并通过数值试验得出了很好的数值结果，从而使时间步长不再受显式的控制，实现了无条件稳定。
【总页数】5页(P80-84)
【作者】向新华;史汉生
【作者单位】解放军理工大学气象学院;解放军理工大学理学院，南京211101【正文语种】中文
【中图分类】O241.82
【相关文献】
1.求解一维抛物型方程的一类半隐格式及组显格式 [J], 姜蕴芝;蔡志权;葛永斌;
2.求解一维抛物型方程的一类半隐格式及组显格式 [J], 姜蕴芝;蔡志权;葛永斌
3.一种时间分数阶扩散方程初边值问题的隐式有限差分格式 [J], 陈春华;卢旋珠
4.求解运动方程的一种不等时间步长的显式数值积分方法 [J], 周正华;王宇欢;刘泉;尹晓涛;杨程
5.一种求解有交界面的椭圆型方程的隐式间断条件差分格式 [J], 徐建军;
因版权原因，仅展示原文概要，查看原文内容请购买。
用向后euler法求三维抛物型方程的adi差分格式
用向后Euler法求解三维抛物型方程的ADI差分格式在数学和科学工程中，三维抛物型方程是一类常见的偏微分方程。
为了求解这种类型的方程，可以使用ADI（Alternating Direction Implicit）差分格式，并结合向后Euler法进行数值计算。
本文将详细介绍如何使用向后Euler法求解三维抛物型方程，并给出相应的ADI差分格式。
1. 三维抛物型方程首先，我们来定义三维抛物型方程。
一个典型的三维抛物型方程可以表示为：∂u/∂t = ∇²u + f(x, y, z, t)其中，u(x, y, z, t)是未知函数，f(x, y, z, t)是已知函数，∇²是拉普拉斯算子。
2. 向后Euler法向后Euler法是一种常用的数值方法，用于求解偏微分方程的时间演化过程。
它基于离散化时间步长，并采用向后差分近似导数。
对于一个一阶时间导数∂u/∂t，向后Euler法将其近似为：(du/dt)_{n+1} ≈ (u_{n+1} - u_n) / Δt其中，u_n和u_{n+1}分别表示时间步n和n+1时刻的解，Δt是时间步长。
3. ADI差分格式ADI（Alternating Direction Implicit）差分格式是一种常用的数值方法，用于求解偏微分方程。
它将空间导数通过交替方向进行离散化，并使用隐式格式来处理。
对于三维抛物型方程，可以使用ADI差分格式来近似求解。
具体步骤如下：•将空间域划分为网格点(x_i, y_j, z_k)，其中i = 0, 1, ..., M_x - 1，j = 0, 1, ..., M_y - 1，k = 0, 1, ..., M_z - 1。
•将时间域划分为时间步n和n+1。
•在每个时间步中，按照以下顺序进行计算：–在 x 方向上进行隐式差分：A_x u^{n+1/2}_{i,j,k} = B_x u^n_{i,j,k} + Δt f^{n+1/2}_{i,j,k}–在 y 方向上进行隐式差分：A_y u^{n+1}_{i,j,k} = B_y u^{n+1/2}_{i,j,k} + Δt f^{n+1}_{i,j,k}–在 z 方向上进行隐式差分：A_z u^{n+3/2}_{i,j,k} = B_z u^{n+1}_{i,j,k} + Δt f^{n+3/2}_{i,j,k}•重复以上步骤直到达到所需的时间步数。
抛物型方程的差分方法
抛物型方程的差分方法抛物型方程是描述物理现象中的薄膜振动、热传导、扩散等过程的方程，具有非常重要的应用价值。
差分方法是一种常用的数值计算方法，用于求解微分方程，对于抛物型方程的数值求解也是非常有效的方法之一、本文将介绍抛物型方程的差分方法，并具体讨论用差分方法求解抛物型方程的一些具体问题。
首先，我们来介绍一下抛物型方程的一般形式。
抛物型方程一般可以表示为：∂u/∂t=α(∂^2u/∂x^2+∂^2u/∂y^2)其中，u(x,y,t)是待求函数，t是时间，x和y是空间变量，α是常数。
这个方程描述的是物理过程中的扩散现象，如热传导过程、溶质的扩散过程等。
差分方法的基本思想是将求解区域离散化为一个个网格点，然后在每个网格点处用近似的方式来计算待求函数的值。
差分方法的求解步骤主要包括以下几个方面：1.选择适当的网格和步长。
在求解抛物型方程时，需要确定空间变量x和y所在的网格点以及步长，同时也需要确定时间变量t所在的网格点和步长。
通常，我们会选择均匀网格，步长选择合适的值。
2.建立差分格式。
差分格式是差分方法的核心部分，它包括对方程进行近似处理和离散化。
对于抛物型方程，常用的差分格式有显式差分格式和隐式差分格式等。
其中，显式差分格式的计算速度快，但是有一定的稳定性限制，而隐式差分格式的稳定性较好，但是计算量较大。
因此，在具体问题中需要根据实际情况选择适当的差分格式。
3.编写计算程序。
在建立差分格式后，需要编写计算代码来求解离散方程。
具体编写的过程包括定义初始条件、建立迭代计算过程、以及计算结果的输出等。
4.计算结果的验证与分析。
求解方程后，需要对计算结果进行验证和分析，主要包括对数值解和解析解的比较、对误差的估计和控制等。
在具体求解抛物型方程时，还会遇到一些问题，例如边界条件的处理、稳定性和收敛性的分析等。
下面将对其中一些问题进行详细讨论。
1.边界条件的处理。
边界条件对差分格式的求解结果有着重要的影响，常见的边界条件包括固定端（Dirichlet）边界条件和自由端（Neumann）边界条件等。
二维抛物型方程的交替方向隐格式
二维抛物型方程的交替方向隐格式二维抛物型方程的交替方向隐格式是一种数值解法，用于求解二维抛物型方程的数值解。
这种解法将二维问题分解为两个一维问题，并采用隐式差分方法来求解。
具体来说，二维抛物型方程可以表示为：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值，我们可以逐步逼近方程的数值解。
这种解法具有二阶精度和稳定性，可以应用于多种二维抛物型方程的问题。
四阶抛物方程一类新的并行交替分段隐格式
四阶抛物方程一类新的并行交替分段隐格式解决四阶抛物方程的并行交替分段隐格式是一种新的数值方法。
被广泛用于研究计算物理，热力学，传热热流等传热流动的复杂问题，是计算科学研究的必要基础。
一、四阶抛物方程的定义四阶抛物方程是一种特殊的常微分方程，它由具有四个平方型、三次项、二次项和定数项的多项式组成，用解析函数求解，有其特殊的几何意义，是研究抛物线在真实世界中几何变化规律的有力工具。
二、并行交替分段隐格式并行交替分段隐格式是一种并行求解抛物方程的方法，与传统的单池迭代法相比，充分利用了计算机的并行计算能力，从而提高了求解效率。
该算法按数值积分的原理，将抛物方程的空间区域划分成多段，每段隐格式均把抛物方程的复杂性分解，实现了对抛物方程的多段重叠解算。
同时，并行交替分段隐格式还使得数值解算过程更加简易，更加精确有效。
三、新颖性1、并行交替分段隐格式利用了技术上的突破，将一个抛物方程划分多段，通过分解求解，从而提高了求解的速度和准确度；2、并行交替分段隐格式在处理抛物方程的空间和时间变化时具有稳定性、快速性及准确性，是解决复杂问题及计算物理、热力学、传热、热流动等传热流动问题的新方法；3、并行交替分段隐格式将原有的抛物方程分为更多段，计算方向不再受求解衔接点的限制，有效解决了抛物方程在均匀网格下的求解困难；4、并行交替分段隐格式在并行计算中有着独特的吸引力，优化了计算的过程，使求解过程更加快捷高效；5、这一算法有时会被称为“交替循环”，它将包含有多个隐格式段的数值求解过程进行方向和分类的循环，从而大大简化了求解过程，提高了求解效率。
四、应用前景随着计算机技术的不断发展和技术成熟，并行交替分段隐格式在危机场景下应用越来越多。
它不仅可以解决抛物方程，还能高效解决求解多维空间与时间变化物理问题。
这种方法同样也能应用在计算风流场、流体动力学方面，也在气动学、建模领域有着广泛的应用。
未来，并行交替分段隐格式将在多学科领域发挥大的作用，促进各项科学研究，为我们创造更加美好的未来。
求解一维抛物型方程的一类半隐格式及组显格式
求解一维抛物型方程的一类半隐格式及组显格式一维抛物型方程是数学中比较常见的一类方程，用来描述各种函数的性质，形式如下：$$ y=ax^2+bx+c $$半隐格式与组显格式是求解一维抛物型方程的两种经典技术，都能应用于解决一维抛物型方程中实数根的问题。
下面就分别介绍这两种技术：1. 半隐格式：半隐格式也称半陈述格式，是一类快速求解一维抛物型方程根的数值求解方法，可以给出一维抛物型方程的根，不必经过复杂的运算。
把一维抛物型方程$ax^2+bx+c=0$改为如下半隐格式：$$ \frac{x^2+\frac{b}{a}x+\frac{c}{a}=0}{a} $$显然，上式的两个根都是$\frac{-b \pm \sqrt{b^2-4ac}}{2a}$，需要注意的是，对于无解的情况，半隐格式也能给出无解的结果。
2. 组显格式：组显格式也称组文格式，是一类快速求解一维抛物型方程实数根的数值求解方法，可以给出一维抛物型方程根，不必经过复杂的运算。
将一维抛物型方程$ax^2+bx+c=0$改为如下组显格式：$$\left\{\begin{aligned}x^2+\frac{b}{a}x+\frac{c}{a}=0 \\x+\frac{c}{a}=0\end{aligned}\right.$$组显格式的两个根分别是$\frac{-b}{2a}$和$\frac{-c}{a}$。
由组显格式可以发现，由于组显格式具有更加直观的单实根示意，因此可以更加快速精确地获得一维抛物型方程的实数根。
总结来说，半隐格式和组显格式都是求解一维抛物型方程的常用技术，它们具有计算迅速、简单直观等优点，半隐格式能求出所有一维抛物型方程的根，而组显格式能够更迅速精确地求出实数根。
它们也常常被应用于应用于微积分、几何学等图形函数的计算。
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方向更新，或者反之）和隐格式（例如，使用向后差分格式近似二阶导数项），将二维抛物型方程中的空间导数进行近似。
通过交替迭代求解这组离散代数方程，我们可以得到二维抛物型方程的数值解。
当离散网格点的数量足够多时，数值解将趋近于方程的解。
解抛物型方程的一族六点隐式差分格式
解抛物型方程的一族六点隐式差分格式詹涌强【期刊名称】《安徽大学学报（自然科学版）》【年(卷),期】2012(036)004【摘要】Proposed in the paper was a class of two-level implicit difference schemes for solving one-dimension parabolic equation. The truncation error of the scheme was O(Τ2 + h4 ) . By Fourier method, the difference scheme was proved to be unconditionally stable if 1/2 ≤ θ ≤ 1, and the stability condition was 0 < r≤ 1/6(1 - 2θ) while 0 ≤ θ < 1/2. The numerical experiment showed the difference scheme was effective and theoretical analysis of them coincides with practical calculation of them.%提出了求解一维抛物型方程的一族两层六点隐式格式.格式的截断误差为O(τ2+h4).利用Fourier 方法证明了差分格式当1/2≤θ≤1时,格式绝对稳定；当0≤θ＜1/2时,只有r≤1/6(1-2θ),格式才是稳定的.数值试验表明,该族格式是有效的,且理论分析与实际计算相吻合.【总页数】4页(P26-29)【作者】詹涌强【作者单位】华南理工大学广州学院基础部数学教研室,广东广州 510800【正文语种】中文【中图分类】O241.82【相关文献】1.解抛物型方程的八点隐式差分格式 [J], 周敏;高学军;董超2.抛物型方程的一族高精度恒稳定的隐式差分格式 [J], 马明书;王同科;申培萍3.抛物型方程一族高精度隐式差分格式AGEI方法 [J], 于战华;金承日4.解高阶抛物型方程的一族隐式差分格式 [J], 曾文平5.解抛物型方程的一族高精度隐式差分格式 [J], 詹涌强;张传林因版权原因，仅展示原文概要，查看原文内容请购买。
抛物型微分方程
数学与计算科学学院实验报告
实验项目名称抛物型方程数值解
所属课程名称微分方程数值解法
实验类型验证
实验日期
班级信计0902
学号
姓名
成绩
{
double uu;
uu=exp(pow(pi,2)*(-t))*sin(pi*x);
return uu;
}
2，结果分析：
根据上面的实验原理，可以得出：在每一时间层，需要求解的隐式差分方程形成了一个线性代数方程组，它的系数矩阵是三对角形矩阵，即仅在主对角线及其相邻二条对角线上有非零元素，故可以使用三角分解法。
使用Crank_Nicolson隐式法，每时间层包含较多的计算工作量，但是其优点为其稳定性要求对步长比的限制大为放宽，而这正是我们所期望的。
【实验结论】（结果）
下面仅显示X轴上取x=0.5处时间轴T的一列间断点出的函数值：
附录1：源程序
附录2：实验报告填写说明。
抛物型方程的O(h 4)精度新交替分段显隐格式
（ ．Ｓｈｏ ｏ ￣ ｌ ｔｓ ｎ ｙｔ ｉ ｃｓ ｈｎｏｇＵ ｉ ｒｔ， ｉ ｎ２００ ，Ｓａｄｎ ，Ｃ ｉ ； １ ｃｏｌ ｆ ｌｎ ｉ ｄＳｓ ｍ Ｓ ｅ ｅ，Ｓａｄｎ ｎ ｅ ｉ Ｊ ａ １０ ｈｎｏｇ ｈｎ Ｍａ ｅ ａｃ ａ ｅ ｃ ｎ ｖ ｓｙ ｎ ５ ａ
ｗ ｉｈ Ｗａ ｒｖ ｄ ｔ ｅ ｈｇ ｌ ｅ ｕ ａ ｙ ｎ ｍｅ ｉａ ｅ ｐ ｒ ｎ ｓ ｈ ｃ ｓ ｐ ｏ ｅ ｂ ｉ ｙ ａ ｃ ＿￣ ｂ ｕ ｒ ｌ ｘ ｅｉ ｏ ｈ ｒｅ ｃ ｍｅ ｔ．
Ｋｅ ｙｗ
：８ １ ａｔ ， ｉ ｓ ｇ ｅ ｔｓｈ ｍｅ ｏ ｖ ｃｉｎ ｄｆ ｓ ｎ ｐ ｂｅ １￣ ｌ ｅ ￣ ， ｅ￣ ｎ ｃ ｅ ；ｃ ｎ ｅ ｔ ｉｕ ｉ ｒｌ ｍ；ｐ ｒｌ ｌ ｌｏｉ ｍ １ ｒ ￣ ｎ ｏ ｏ ｏ ａａｌ ｇ ｒ ｈ ｅａ ｔ
ｃ ｍ ｎ ｅ ａ ｌｌ ｍ ｕ ｔ ｎ ｗ ｉ ａ ｅ ｅｕｃｎｉｏ ｌ ａｌ． ｈ ｆ ｕ ｃ ｔ ｉ Ｏ ｈ ） ｓ ｅ ｅ ａ ｅ ｓｄｉ ｐｒｌ ｏ ｐｔ ｏ ， ｈ ｈｗｓ ｒｖｄｔ ｂ ｎｏｄ ｉ ａｙｓｂｅ Ｔ ｅｓａＭ ａｃｒ ｙ ａ （ ， ｈ ｃ ｂｕ ｎ ａｅ ｃ ａ ｉ ｃ ｐ ｏ ｏ ｔｎｌ ｔ ｐｉ ｃ ａ ｒ ｅ ｓ
维普资讯
第４卷 ２
Ｖ０ ． ２ １４
第 ｌ 期 ２
Ｎ ．２ ｏ １
山
东
大
学
学
报
（ 理
学
版）
２０ 年 ｌ Ｏ７ ２月
Ｄｅ ｃ．２ ｒ ００７
Ｊｕｎｌ ｆ ｈｎ ｏｇＵ ｉ ｒ ｔ（ ａ ｒ ｃｅｃ） ｏｒａ ｏ ａｄｎ ｎｖ ｓｙ Ｎ ｔ ａ Ｓｉ ｅ Ｓ ｅｉ ｕｌ ｎ
基于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.
偏微分方程数值解法(抛物型方程差分法)2
显格式 隐格式
1
[u
k j
1
u
k j
]
a2 h2
x2u
k j
1
[u
k j
1
u
k j
]
a2 h2
u2 k1 xj
C-N格式
1
[u
k j
1
u
k j
]
a2 2h2
2 x
(u
k j
1
u
k j
)
15/17
数值计算实验 显格式: input T:=1 error = 7.9443e-006 k = 200
C—N格式 input T:=1 error = 2.6227e-006 k = 50
12/17
C-N 格式矩阵形式
[(1 ra2 )I ra2 C ]uk1
[(1
ra 2
2 )I
ra 2
C ]uk
(fk
f k1 )
2
2
H [(1 ra2 )I ra2 C ]1[(1 ra2 )I ra2 C ]
2
2
特征值
j
(1 (1
ra2 ) ra2 )
ra 2 ra 2
cos( cos(
•
8、业余生活要有意义，不要越轨。20 20年12 月10日 星期四 5时54 分52秒0 5:54:52 10 December 2020
•
9、一个人即使已登上顶峰，也仍要自 强不息 。上午 5时54 分52秒 上午5时 54分05 :54:522 0.12.10
• 10、你要做多大的事情，就该承受多大的压力。12/10/
( j 1,2,, n)
|
数学中的抛物型方程
数学中的抛物型方程抛物型方程（parabolic equation）是数学中一类重要的偏微分方程，它在物理学、工程学和社会科学等领域中具有广泛的应用。
本文将从抛物型方程的定义、特征和解法等方面进行论述，以帮助读者更好地理解和应用抛物型方程。
一、抛物型方程的定义在数学中，抛物型方程是一类二维或三维偏微分方程，其形式可以表示为：∂u/∂t = a∇²u + bu + c其中，∂u/∂t 表示函数 u 对时间 t 的偏导数，∇²u 表示函数 u 对空间坐标的拉普拉斯算子，a、b、c 是常数。
抛物型方程通常描述了某一物理现象随时间变化的规律，比如热传导、扩散等。
通过解抛物型方程，我们可以预测和分析这些物理现象。
二、抛物型方程的特征1. 热传导方程抛物型方程在热传导方程中的应用是最常见的。
热传导方程描述了物体内部温度随时间和空间的变化情况。
在一维情况下，热传导方程具有以下形式：∂u/∂t = α∂²u/∂x²其中，u(x, t) 表示在时刻 t 位置为 x 的温度，α 是热扩散系数。
2. 扩散方程抛物型方程在扩散方程中的应用也是非常重要的。
扩散方程描述了物质在浓度梯度驱动下的扩散过程。
在一维情况下，扩散方程具有以下形式：∂u/∂t = D∂²u/∂x²其中，u(x, t) 表示在时刻 t 位置为 x 的物质浓度，D 是扩散系数。
三、抛物型方程的解法对于抛物型方程，我们通常采用偏微分方程的求解方法，如分离变量法、格林函数法等。
1. 分离变量法分离变量法是一种常用的求解抛物型方程的方法。
它的基本思想是将多元函数分解为几个一元函数的乘积，并利用分离后的一元函数满足各自的方程来求解。
以热传导方程为例，我们可以将其分离变量为时间部分和空间部分：u(x, t) = X(x)T(t)代入原方程，得到两个方程：X''(x)T(t)/X(x) = T'(t)/T(t) = -λ²其中，λ² 是常数。
有限差分法求解抛物型方程说明
有限差分法求解抛物型方程偏微分方程只是在一些特殊情况下，才能求得定解问题解的解析式，对比较复杂的问题要找到解的解析表达式是困难的，因此需采用数值方法来求解.有限差分法是一种发展较早且比较成熟的数值求解方法，只适用于几何形状规则的结构化网格.它在微分方程中用差商代替偏导数，得到相应的差分方程，通过解差分方程得到微分方程解的近似值.本章主要介绍有限差分法的基本思想，并给出一些具体的数值实例.§1 差分方法的基本思想有限差分法把偏微分方程的求解区域划分为有限个网格节点组成的网格，主要采用Taylor 级数展开等方法，在每个网格节点上用有限差分近似公式代替方程中的导数，从而建立以网格节点上的函数值为未知数的代数方程组.有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式.从差分的空间形式来考虑，可分为中心格式和逆风格式.考虑时间因子的影响，差分格式还可以分为显格式、隐格式和显隐交替格式等.目前常见的差分格式，主要是上述几种格式的组合，不同的组合构成不同的差分格式.泰勒级数展开法对有限差分格式的分类和公式的建立起着十分重要的作用.下面采用泰勒展开式导出一个自变量系统的若干有限差分表达式.首先考虑单变量函数()u x ，如图1把区域x 离散为一批结点，记0()(), =0,1,2,i i u x u x ih u i =+=图1 单变量函数离散化函数()u x 在点i x 处的泰勒展开式为23()()()()()2!3!i i i i i u x u x u x h u x u x h h h ''''''+=++++ （1）或23()()()()()2!3!i i i i i u x u x u x h u x u x h h h ''''''-=-+-+ （2）式(1)和(2)重新整理可得2()()()()()2!3!i i i i i u x h u x u x u x u x h h h '''''+-'=---（3）和2()()()()()2!3!i i i i i u x u x h u x u x u x h h h '''''--'=+++（4）于是给出在点i x 处函数u 的一阶导数的两个近似公式1()()()i i i ii u x h u x u u u x h h ++--'≈= （5）1()()()i i i i i u x u x h u u u x h h----'≈= （6）因为级数被截断，这两个近似公式肯定要产生误差，此误差与h 同阶，形式分别为()(), ,2()(), .2i i i i i i hE u O h x x h hE u O h x h x ξξξξ''=-=≤≤+''==-≤≤ 若把式(3)和(4)相加并求()i u x '，可得11()()()22i i i i i u x h u x h u u u x h h+-+---'≈= （7）其截断误差与2h 同阶，形式为22()(), ,6i i i h E u O h x h x h ξξ''=-=-≤≤+若把式(3)和(4)相减并求()i u x ''，可得1122()2()()2()i i i i i i i u x h u x u x h u u u u x h h +-+-+--+''≈= （8）其截断误差与2h 同阶，其形式为22()(), ,12i i i h E u O h x h x h ξξ''=-=-≤≤+我们可继续用这种方式来推导更复杂的公式，类似的公式还有很多，这里不再一一列举.公式（5）、（6）分别称为一阶向前、向后差分格式，这两种格式具有一阶计算精度，公式（7）、（8）分别称为一阶、二阶中心差分格式，这两种格式具有二阶计算精度.图2 二维区域网格剖分上面的结果可直接推广使用于导出二元函数(,)u x y 的许多有限差分近似公式.如图7.2，把求解区域进行网格剖分，使12(,)(,), ,=0,1,2,i j ij u x y u ih jh u i j ==其中x 方向的网格间距为1,h y 方向的网格间距为2,h 整数i 和j 分别表示函数(,)u x y 沿x 坐标和y 坐标的位置.二元函数(,)u x y 对x 求偏导时y 保持不变，对y 求偏导时x 保持不变，根据向前差分公式（7.5）可以给出在点(,)i j x y 处函数(,)u x y 的一阶偏导数的两个近似公式1,,1(,)i j i j i ju x y u u xh +∂-≈∂ （9）,1,2(,)i j i j i ju x y u u yh +∂-≈∂ （10）相类似地，根据二阶中心差分格式（8）可以得到函数(,)u x y 的二阶偏导数的近似公式21,,1,221(,)2i j i j i j i ju x y u u u x h +-∂-+≈∂ （11）2,1,,1222(,)2i j i j i j i j u x y u u u yh+-∂-+≈∂ （12）下面我们推导函数(,)u x y 的二阶混合偏导数2ux y∂∂∂在(,)i j x y 的有限差分表达式.根据一阶中心差分格式（7），112111,11,11,11,122121221,11,1(,)(,)(,)1()21 ()()222 i j i j i j i j i j i j i j i j i j i u x y u x y u x y O h x y h y y u u u u O h O h h h h u u u +-+++--+--+++-∂∂∂⎡⎤⎡⎤∂=-+⎢⎥⎢⎥∂∂∂∂⎣⎦⎣⎦--⎡⎤=-++⎢⎥⎣⎦--≈1,11,1124j i j u h h -+--+二维有限差分近似可以直接推广到三维空间或三维空间加一维时间的情形.定义1 当步长趋于零时，差分方程的截断误差趋于零，则称差分格式与微分方程是相容的.定义2 当步长趋于零时，差分方程的解收敛于微分方程的解，则称差分格式是收敛的. 定义3 当差分方程的解由于舍入误差的影响，所产生的偏差可以得到控制时，则称差分格式是稳定的.§2 抛物型方程的有限的差分法为了说明如何使用有限差分法来求解偏微分方程，本节我们给出以下几个数值实例.算例1 考虑一维非齐次热传导方程的初边值问题：2212(,), 01,01,(,0)(), 01,(0,)(), (1,)(), 0 1.u ua f x t x t t x u x q x x u t g t u t g t t ⎧∂∂=+<<<≤⎪∂∂⎪⎪=≤≤⎨⎪==<≤⎪⎪⎩(7.13),其中2,a =函数11(,)[cos()2sin()],22xf x t e t t =--+-初始条件1()sin,2xq x e =左、右边界条件分别为11()sin(),2g t t =-21()sin()2g t e t =-.该定解问题的解析解为1(,)sin(),(,)[0,1][0,1].2xu x t e t x t =-∈⨯将求解区域{(,)|,0}x t a x b t T Ω=≤≤≤≤进行网格剖分，[,]a b 作m 等分，[0,]T 作n 等分，记,,b a Th m nτ-==则 ,0,,0i k x a ih i M t k k n τ=+≤≤=≤≤对该问题建立如下向前差分格式：11122, 11, 11,k kk k k k i i i i i i u u u u u a f i m k n hτ+-+--+=+≤≤-≤≤-（14） (,0)(),1,i i u x q x i m =≤≤ （15） 12(,)(), (,)(),1.k k k k u a t g t u b t g t k n ==≤≤ （16）令2r ah τ=，差分格式（7.14）整理得111(12), 11, 1 1.k k k k k i i i i i u ru r u ru f i m k n τ+-+=+-++≤≤-≤≤- （17）显然时间在1k t +上的每个逼近值可独立地由k t 层上的值求出。
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 取值。
-解抛物型方程的一个新的六点隐格式
圳0<r≤ 1 2
3 数值例子
对初边值问题：
坠u 坠t
=
坠2u 坠x2
,0≤x≤1.0,t>0
u(x,0)=sinx,0≤x≤1.0
（8）
u(0,t)=0,u(1,t)=exp(-t)sin(1,0),t≥0
已知问题（8）的精确解为 u(x,t)=e-1sinx，我们利用格式（7）求解数值解，并和精确解进行比较。 表 1
(4)
解(4)式，并整理得
η1=(3r-1)η5
η2=rη5
(5)
η3=rη5
η4=-η5
n
n
取 η5 =1，将 ηi(i=1，2，3，4）代入（3），舍去截断误差，用 uj 表示 u j 的近似值,得差分格式为:
n
n
n
n+1
n+1
n
η1uj+1 +η2uj-1 +η3uj+2 +η4uj +η5uj+1 -(η1+η2+η3+η4+η5)uj
豆丁推荐 ↓
精品文档
曹俊英 1 ， 王自强 1 ， 张大凯 2
（1、贵州民族学院数学与计算机科学学院，贵州 贵阳 550025；2、贵州大学数学系，贵州 贵阳 550025）
摘 要： 用待定系数法对一维抛物型方程构造了一个两层六点差分隐格式,使得精度达到 O (τ2+h3)， 稳定性条件为 0<r ≤1/2。
（责编： 谭本龙 责校： 张永光)
·48·
当 模 型 （1）的 解 充 分 光 滑 时 ，有 如 下 关 系 式 成 立 ：
坠m+nu 坠xm坠tn
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
2. 实验内容：
用 -方法导出的隐式格式求如下初边值问题
u u 0 , 0 x 1, 0 t 1 2 t x u ( x , 0 ) sin x , 0 x 1 , u ( 0 , t ) 0 , u (1 , t ) 0 , 0 t 1 ,
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0001 -0.0053 0.3102
0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0009 -0.0706
U
n j 1
2U
n j 2
U
n j 1
( x)
j
( x j ), U
2
0, U
n J
0
以r 为
t /( x )
表示网比。由于本是实验中
( x ) sin x
, 所以上述差分格式简化
rU U
0 j
n 1 j 1
(1 2 r )U
运行结果如下： cita=0 M=10 J=10 x = 1.0e+010 * Columns 1 through 9 0 0.0099 Column 10 0 -2.4999 3.6685 -3.2575 2.0465 -0.9383 0.3102 -0.0706
u = 1.0e+010 * Columns 1 through 9
5
0.0000 0.0000 0 0.0000 0 0.0000 0 0.0000 0 0.0000 0 0.0000 0 -0.0000 0 0.0000 0 -0.0000 0 -0.0001 0 0.0099
0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0001 -0.0023 0.0754 -2.4999
淮
海
工
学
院
实 验 报 告 书
课程名称： 实验名称： 班 姓 日 级 名： 期：
微分方程数值解法 抛物型方程隐格式的程序实现 信息 101
学号： 地点 成绩：
数学实验室
指导教师：
数 理 科 学 系
1
1. 实验目的：
按照实验内容对 -方法导出的几种隐式差分格式编程实现求出数值解； 估计 数值解的误差界；大致比较几种不同格式的优劣；结合格式的相容性、稳定性和 收敛性条件简单分析计算结果。
间步长为 x 线：
x x
j
1/ J
和时间步长为 t
，其中 J ，M 为正整数，用两族平行直
j x , j 1 , 2 , ..., J
平行于 t 轴作竖直线和 平行于 x 轴作水平线,
t t n n t , n 1 , 2 , ..., M
将矩形区域 G 分割成矩形网格，得网格剖分图。 第二步: 差分法的目的是：求初边值问题（2）的解 u ( x , t ) 在节点 ( x j , t n ) 的近 似值 U
0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0012 -0.0505 2.0465
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0004 0.0200 -0.9383
Columns 10 through 11 0.0000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
h = 1.0e+009 * 0.0000 0.0007 -0.7543 1.0667 -0.8874 0.5047 -0.2001 0.0534 -0.0087
0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0001 0.0031 -0.1067 3.6685
0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0001 -0.0024 0.0887 -3.2575
j
( u ) [(
1 2
) tu
xxt
1 8
1 12
( x) u
2
n xxxx
1 2
]
1 2
j
[
1 24
( t ) u ttt
2
( t) u
2
n xxtt
]
j
O (( t )
3
( x) )
4
利用 Fourier 方法判断差分格式的稳定性条件是： 当0
1 2
，当且仅当 r
1 2
(1 2 )
1
时稳定，当
1 2
1 ，对所有的 r
格式
均稳定。 利用最大值原理可以证明差分格式的收敛条件是： r
1 2 (1 )
1
。
第三步: 有限差分方程的解法。根据问题的规模和计算机的容量和速度，选 取适当的解法。这是数值计算的关键一步，有限差分法解方程的主要计算都集中 在这里。 由于本实验问题的有限差分方程为 -方法导出的隐式差分格式，可以用 Thomas 方法求解。对每个时间层 n 1 需要求解的方程为：
n 0
0, U
n J
0
。
第四步: 根据前面的分析，编制程序，上机计算，求出数值解，必要时画出 解的几何图形，分析观察解的性质。 为了方便比较可以用傅立叶级数方法求出定解问题的解析解为：
u ( x, t) e
t
2
sin x
。
另外，上述差分程组为三对角方程组，可用 Thomas 法（追赶法）求解。 若有三对角方程组：
n j
( j 1, 2 , , J 1; n 1, 2 ,..., M ) 。为此，需要构造逼近微分方程定解问
题的差分格式。采用 -方法，构造如下差分格式
U
n 1 j
U
n jt U0 j NhomakorabeaU
n 1 j 1
2U
n 1 j 2
U
n 1 j 1
( x)
n 0
(1 )
6
>>
cita=0.5 M=10 J=10 x= 1.0e-003 * Columns 1 through 9 0 -0.0037 Column 10 0 -0.9215 0.2009 0.2907 0.0655 -0.0568 -0.0596 -0.0252
7
u= Columns 1 through 9 0.3090 0.3090 0 0.1041 0 0.0305 0 0.0063 0 0.0017 0 0.0006 0 -0.0000 0 0.0001 0 -0.0000 0 0.0000 0 -0.0000 0.5878 0.1335 0.0163 0.0139 -0.0029 0.0037 -0.0023 0.0019 -0.0014 0.0011 -0.0009 0.8090 0.2338 0.0425 0.0174 0.0012 0.0017 -0.0002 0.0001 0.0001 -0.0002 0.0002 0.9511 0.2981 0.0664 0.0183 0.0046 0.0005 0.0009 -0.0005 0.0005 -0.0004 0.0003 1.0000 0.3250 0.0813 0.0185 0.0061 0.0004 0.0009 -0.0003 0.0002 -0.0001 0.0001 0.9511 0.3148 0.0846 0.0178 0.0059 0.0008 0.0005 0.0000 0.0000 0.0001 -0.0001 0.8090 0.2704 0.0760 0.0156 0.0049 0.0010 0.0002 0.0002 -0.0001 0.0001 -0.0001 0.5878 0.1975 0.0571 0.0117 0.0034 0.0009 0.0000 0.0002 -0.0001 0.0001 -0.0000
，
3
xn fn x i f i e i x i 1 , i n 1 , n 2 , ,1
其中 x 1 ,
x 2 , ..., x n
为所求三对角方程组的解。
四． 实验数据记录及分析（或程序及运行结果） ：
程序如下：
clear cita=input('cita='); M=input('M='); J=input('J='); u=zeros(J+1,M+1); d=zeros(1,J); e=zeros(1,J-1); f=zeros(1,J); x=zeros(1,J); uu=zeros(1,J); for n=1:M u(n,1)=0; u(n,M)=0; end for j=1:J u(1,j)=sin(j*pi/J); end tt=1/M; xx=1/J; r=tt/(xx*xx); for n=1:M for j=2:J-1; d(1,j)=(1-cita)*r*u(n,j+1)+(1-2*(1-cita)*r)*u(n,j)+(1-cita)*r*u(n,j-1) ; a=cita*r; b=1+2*cita*r; c=cita*r; end e(1,2)=c/b; f(1,2)=d(1,2)/b; for i=3:J-2 e(1,i)=c/(b-e(1,i-1)*a); end for i=3:J-1