二维抛物线方程数值解法(ADI隐式交替法)方法
高中数学学习中的抛物线与双曲线方程求解方法
高中数学学习中的抛物线与双曲线方程求解方法在高中数学学习中，抛物线与双曲线是重要的二次函数的图像形式。
学生们需要掌握求解抛物线和双曲线方程的方法，以便能够准确地描述并解决与这些图形相关的问题。
本文将介绍高中数学学习中抛物线与双曲线方程求解的方法。
首先，我们来讨论抛物线的方程求解。
一般来说，抛物线的方程通常是二次函数的形式：y = ax² + bx + c。
在求解抛物线方程时，我们通常要考虑以下几种情况：一种情况是已知抛物线上的三个点，我们需要确定抛物线的方程。
对于已知三个点（x₁，y₁）、（x₂，y₂）、（x₃，y₃），我们可以建立三个方程：(1) y₁ = ax₁² + bx₁ + c(2) y₂ = ax₂² + bx₂ + c(3) y₃ = ax₃² + bx₃ + c通过解这个方程组，我们可以找到抛物线的方程。
另一种常见情况是已知抛物线的顶点和一点，需要确定抛物线的方程。
对于已知顶点（h，k）和一点（x₁，y₁），我们可以通过将这两个点代入抛物线的一般方程，得到下面的方程：(1) y₁ = a(x₁ - h)² + k通过解这个方程，我们可以得到抛物线的方程。
在实际问题中，我们常常需要求解与抛物线相关的问题。
例如，给定一个抛物线，我们需要找到它的焦点和准线。
对于抛物线方程 y = ax² + bx + c，我们可以通过求解以下方程得到焦点（p，q）和准线的方程：(1) p = -b / (2a)(2) q = c - (b² - 1) / (4a)通过求解这两个方程，我们可以找到焦点和准线的方程。
接下来，我们转到双曲线的方程求解。
与抛物线类似，双曲线的方程也是二次函数的形式：y = a/x。
在求解双曲线方程时，我们同样需要考虑不同的情况。
一种情况是已知双曲线上的两个点，我们需要确定双曲线的方程。
对于已知两个点（x₁，y₁）和（x₂，y₂），我们可以建立以下方程：(1) y₁ = a/x₁(2) y₂ = a/x₂通过解这个方程组，我们可以找到双曲线的方程。
二阶二维抛物型方程混合问题的样条解法
?
本 文 旨在 解 决 这一 问题
:
郡 将〔3 〕 的方 法 推广 到解 二 阶二 维抛 物型 方 程 的棍
设 二 阶 二 维 抛物 型 方 程 为
L ”二
+
刁U 刁f
一
以劣
,
梦 t)
,
月 2 况 一 。( 二
,
,
,
t )
+ `
月仍
2
哥
(
,
,
(劣
,
“` )
,
会
R
z
d (劣 g t )
,
石 爹
, ,
,
g
,
近 似 地代 替满 足
1)
、
( 1
.
) 的 精 确解 2
为使 叙 述 方 便
4〕 的记 本文沿用〔
o
(o
,
o
,
0 )
.
(o
,
j
,
o )
(o
,
n
,
o )
,
矛
`
l
l 奋
了
夕
l
( i
,
0
,
0)
才一
图 当
t:
二
(
n
。
。 晚 )
0 :
)
又丫丫只 又义只只 丫又火只 火火又又
图
2
n
.
0
)
1
0
时
,
在 正 方 形 口A B C 中 ( 图 2 )
,
由( 1
, , ,
,
2 )
”
可知
` ,
厂 “ { “
二维双曲型方程的隐格式解法
本科毕业论文（设计）题目二维双曲线型方程的交替方向隐格式解法院（系）数学系专业数学及应用数学学生姓名周玲玲学号 10022156指导教师陈淼超职称讲师论文字数 9500完成日期: 2014年6月8日巢湖学院本科毕业论文(设计)诚信承诺书本人郑重声明：所呈交的本科毕业论文(设计)，是本人在导师的指导下，独立进行研究工作所取得的成果。
除文中已经注明引用的内容外，本论文不含任何其他个人或集体已经发表或撰写过的作品成果。
对本文的研究做出重要贡献的个人和集体，均已在文中以明确方式标明。
本人完全意识到本声明的法律结果由本人承担。
本人签名：日期：巢湖学院本科毕业论文 (设计)使用授权说明本人完全了解巢湖学院有关收集、保留和使用毕业论文 (设计)的规定，即：本科生在校期间进行毕业论文(设计)工作的知识产权单位属巢湖学院。
学校根据需要，有权保留并向国家有关部门或机构送交论文的复印件和电子版，允许毕业论文 (设计)被查阅和借阅；学校可以将毕业论文(设计)的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存、汇编毕业，并且本人电子文档和纸质论文的内容相一致。
保密的毕业论文(设计)在解密后遵守此规定。
本人签名：日期：导师签名：日期：二维双曲型方程的交替方向隐格式解法摘要在解决偏微分方程中二维双曲线型方程的问题求解初边值问题时，显示格式的稳定性是有条件的，并且多维的稳定性条件更严，为得到稳定性好的格式，隐式格式受到重视。
用隐式格式求解二维问题得到的线性方程组其系数矩阵为宽带状，因此求解不甚便利，采用交替方向隐式（ADI ）格式可以避免此问题。
本文以基本概念和基本方法为主，同时结合算例实现算法。
第一部分介绍二维双曲线型方程的交替方向隐格式解法的基本概念，引入本文的研究对象——二维波动方程： ，x R ∈，],0(T t ∈第二部分介绍上述方程的二维双曲线型方程的交替方向隐格式及这种格式的存在性、收敛性及稳定性。
六点对称法,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)τ取值。
二阶微分方程数值求解
二阶微分方程数值求解
要数值求解二阶微分方程，首先需要将其转化为一个一阶微分方程组。
假设待求解的二阶微分方程为：
y''(x) = f(x, y(x), y'(x))
将其转化为一阶微分方程组：
z(x) = y'(x)
z'(x) = f(x, y(x), z(x))
然后，可以选择数值求解方法，如欧拉方法、改进的欧拉方法、四阶龙格-库塔方法等等，对这个一阶微分方程组进行数值求解。
以欧拉方法为例，假设已知初始条件 y(x0) = y0，z(x0) = z0，
选择步长 h。
则可以按照以下步骤进行数值求解：
1. 初始化步数 n = 0，设置初始条件 y(x0) = y0，z(x0) = z0。
2. 计算下一步的值：y(x + h) = y(x) + h * z(x)，z(x + h) = z(x) +
h * f(x, y(x), z(x))。
3. 将 x 增加 h，即 x = x + h。
4. 将步数 n 增加 1，即 n = n + 1。
5. 重复步骤 2-4，直到达到目标位置的 x 值（如终点 x 结束条
件 x >= x_end）。
需要注意的是，数值求解方法的精度和稳定性都会受到步长的影响，过大的步长可能导致数值不稳定，过小的步长可能导致
计算量过大。
因此，选择合适的步长是很重要的。
值得一提的是，当二阶微分方程为边值问题时，可以采用有限差分法、有限元法等数值方法进行数值求解。
这些方法会更为复杂，并涉及到边界条件的处理。
二维抛物型方程的交替方向隐格式
二维抛物型方程的交替方向隐格式二维抛物型方程的交替方向隐格式是一种数值解法，用于求解二维抛物型方程的数值解。
这种解法将二维问题分解为两个一维问题，并采用隐式差分方法来求解。
具体来说，二维抛物型方程可以表示为：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.判断抛物线的开口方向。
可以通过方程的二次项系数的符号来判断抛物线的开口方向。
当二次项系数大于零时，抛物线开口向上；当二次项系数小于零时，抛物线开口向下。
三、求抛物线与坐标轴交点通过解方程来求解抛物线与坐标轴的交点，这是很常见的题型。
有两种常用的方法。
1.因式分解法。
将抛物线的方程进行因式分解后，可以得到解析解或根的个数。
进一步，通过观察与分析，可以得出与坐标轴交点的具体坐标。
2.二次函数求根公式。
通过应用二次函数求根公式，可以得到抛物线与坐标轴交点的解析解。
需要注意的是，二次函数求根公式只适用于已经化为标准形式的抛物线。
四、求抛物线的切线与法线求抛物线的切线与法线是一类较难的题型，需要熟练掌握相关的知识与求解方法。
下面将介绍两种常见的方法。
1.切线与法线的斜率法。
通过斜率法可以求得切线与法线的斜率表达式。
具体而言，对于给定的抛物线方程，我们可以通过计算其导数来求得切线或法线的斜率表达式，然后利用该斜率表达式求解切线或法线的方程。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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。