下载文档原格式
本科毕业论文（设计）题目二维双曲线型方程的交替方向隐格式解法院（系）数学系专业数学及应用数学学生姓名周玲玲学号 10022156指导教师陈淼超职称讲师论文字数 9500完成日期: 2014年6月8日巢湖学院本科毕业论文(设计)诚信承诺书本人郑重声明：所呈交的本科毕业论文(设计)，是本人在导师的指导下，独立进行研究工作所取得的成果。
除文中已经注明引用的内容外，本论文不含任何其他个人或集体已经发表或撰写过的作品成果。
对本文的研究做出重要贡献的个人和集体，均已在文中以明确方式标明。
本人完全意识到本声明的法律结果由本人承担。
本人签名：日期：巢湖学院本科毕业论文 (设计)使用授权说明本人完全了解巢湖学院有关收集、保留和使用毕业论文 (设计)的规定，即：本科生在校期间进行毕业论文(设计)工作的知识产权单位属巢湖学院。
学校根据需要，有权保留并向国家有关部门或机构送交论文的复印件和电子版，允许毕业论文 (设计)被查阅和借阅；学校可以将毕业论文(设计)的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存、汇编毕业，并且本人电子文档和纸质论文的内容相一致。
保密的毕业论文(设计)在解密后遵守此规定。
本人签名：日期：导师签名：日期：二维双曲型方程的交替方向隐格式解法摘要在解决偏微分方程中二维双曲线型方程的问题求解初边值问题时，显示格式的稳定性是有条件的，并且多维的稳定性条件更严，为得到稳定性好的格式，隐式格式受到重视。
用隐式格式求解二维问题得到的线性方程组其系数矩阵为宽带状，因此求解不甚便利，采用交替方向隐式（ADI ）格式可以避免此问题。
本文以基本概念和基本方法为主，同时结合算例实现算法。
第一部分介绍二维双曲线型方程的交替方向隐格式解法的基本概念，引入本文的研究对象——二维波动方程： ，x R ∈，],0(T t ∈第二部分介绍上述方程的二维双曲线型方程的交替方向隐格式及这种格式的存在性、收敛性及稳定性。
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值，我们可以逐步逼近方程的数值解。
这种解法具有二阶精度和稳定性，可以应用于多种二维抛物型方程的问题。
关于解二维椭圆型高精度差分方程的交替方向迭代法交替方向迭代法（alternating direction iterative method，简称ADI法）是一种常用于解二维椭圆型高精度差分方程的数值方法。
它的基本思想是将二维问题拆分为两个一维问题，然后依次迭代求解，从而达到解二维问题的目的。
一、AD方法的原理二维椭圆型方程可以表示为：>(A1u(i+1,j)+B1u(i,j+1)-(2A1+B1+C1)u(i,j)+B1u(i,j-1)+A1u(i-1,j))/Δx^2+(D1Δy^2)u(i,j)=f(i,j)，式中i,j分别为空间坐标，f(i,j)为已知的外力场，D1为系数。
对于前一半时间步，我们考虑j方向上的一维问题，可将上述方程改写为：>(A1+B1)u(i+1,j)+(C1-2A1-B1)u(i,j)+A1u(i-1,j)=f1(i,j)，式中f1(i,j)为已知。
对于后一半时间步，我们考虑i方向上的一维问题，可将上述方程改写为：>(A2+B2)u(i,j+1)+(C2-2A2-B2)u(i,j)+A2u(i,j-1)=f2(i,j)，式中f2(i,j)为已知。
通过交替迭代这两个方程，逐步逼近真实解。
二、AD方法的步骤1.选取初始解u0(i,j)，以及参数Δx，Δy和迭代终止准则ε。
2.开始迭代，设定初始误差e0=ε+1，n=0。
3. 使用前一半时间步的方程求解每个i所对应的j方向的一维问题，得到一个临时解ut(i,j)。
4.使用后一半时间步的方程求解每个j所对应的i方向的一维问题，得到一个新的解u(n+1)(i,j)。
5. 计算误差en=，u(n+1)-u(n)，2，n=n+16. 如果en ≤ ε，迭代结束，得到近似解u(n+1)；否则返回第3步。
7.算法结束。
三、AD方法的优势1.AD方法是一种分离的迭代方法，当网格较大时其计算量相对较小。
2.AD方法是一种隐式迭代方法，对稳定性和收敛性有较好的保证。
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。
二维抛物型方程的交替方向隐格式在数学领域中，二维抛物型方程是一类重要的偏微分方程，它们在众多实际问题的数学建模中起着关键作用。
对于这类方程，交替方向隐格式是一种常用且有效的数值求解方法。
本文将详细介绍二维抛物型方程及其交替方向隐格式的原理和应用，希望能为读者提供一份生动、全面且有指导意义的参考材料。
首先，我们来了解什么是二维抛物型方程。
通常，二维抛物型方程可以表示为以下形式：∂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方向更新，或者反之）和隐格式（例如，使用向后差分格式近似二阶导数项），将二维抛物型方程中的空间导数进行近似。
通过交替迭代求解这组离散代数方程，我们可以得到二维抛物型方程的数值解。
当离散网格点的数量足够多时，数值解将趋近于方程的解。
二维抛物线方程数值解法方法二维抛物线方程是数学中常见的偏微分方程之一，表示一个二维区域中的传热、扩散、波动等现象。
求解二维抛物线方程通常需要使用数值解法，其中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方向）的方程组中求解。
通过交替进行，可以逼近二维抛物线方程的解。
解二阶抛物型方程含参数高精度两层差分格式引言在数值计算中，抛物型偏微分方程是一类重要的方程，涉及到热传导、扩散、反应扩散等许多实际问题的数值模拟。
而解二阶抛物型方程含参数的问题则更具有挑战性，建立高精度的差分格式成为研究的重点之一。
本文将详细探讨解二阶抛物型方程含参数的高精度两层差分格式。
二阶抛物型方程含参数先考虑一个边界平滑、参数依赖的二阶抛物型方程：其中，u(x,t)是未知函数，a(x)是参数函数，f(x,t)是外力项。
本文的目标是构建一个高精度的差分格式来求解这个方程。
高精度差分格式的结构由于我们要构建高精度的差分格式，自然而然地思路是采用两层差分格式的结构。
两层差分格式具有更高的精度和更好的稳定性，对于求解复杂的偏微分方程尤为有效。
一维差分格式首先，我们先介绍一维的差分格式。
对于一维问题，我们可以将区域进行离散化，将连续的问题转换为离散的问题。
一维差分格式可以表示为：这个差分格式包含三个部分，分别是时间项、空间项和外力项。
其中，时间项计算了时间导数，空间项使用中心差分方法计算了空间导数，外力项贡献了外力。
多维差分格式对于二维问题，我们可以将其推广为多维差分格式。
多维差分格式的基本思想是利用乘积形式构造更高维的差分格式。
在这里，我们只考虑二维情况：这个差分格式在时间、x方向和y方向都有相应的差分项。
同样，时间项计算了时间导数，空间项使用了中心差分方法计算了空间导数，外力项贡献了外力。
参数依赖的高精度两层差分格式在介绍了一维和多维差分格式后，我们接下来考虑含参数的高精度两层差分格式。
对于含参数的问题，我们需要在空间项中考虑该参数的影响。
高精度的两层差分格式可以表示为：其中，αi,j n是参数依赖的系数。
高精度差分格式的稳定性和收敛性分析在构建高精度差分格式时，我们不仅需要关注其精度，还需要考虑其稳定性和收敛性。
在这一部分，我们将对所构建的高精度差分格式进行稳定性和收敛性的分析。
稳定性稳定性是差分格式是否收敛的重要条件之一。
解二阶抛物型方程含参数高精度两层差分格式
解二阶抛物型偏微分方程是许多领域中非常重要的问题，例如热传导、扩散、化学反应等。
然而，在实际计算中，由于参数的存在，往往需要使用高精度的差分格式来保证数值计算的准确性。
以下是一种含参数的高精度两层差分格式，可以用于解决二阶抛物型偏微分方程：
1. 先验估计
为了使用高精度两层差分格式，我们首先需要进行先验估计，以确定合适的时间步长和空间步长。
具体来说，我们可以使用稳定性分析来确定时间步长和空间步长的上限值。
2. 差分格式
在确定了时间步长和空间步长之后，我们可以开始使用高精度两层差分格式来求解二阶抛物型偏微分方程。
该差分格式通常包括以下几个步骤：
（1）先用向前差分公式求解第一层，得到一个中间解。
（2）再采用Crank-Nicolson格式对第二层进行差分，同时使用前一步得到的中间解进行修正。
（3）最后，将得到的数值解反推回到未知函数的值域中，得到方程的
数值解。
需要注意的是，在使用这种高精度差分格式进行计算时，我们需要使
用高精度的算法来保证计算的准确性。
3. 参数调节
由于实际问题中经常存在参数不确定性的情况，因此，在进行数值计
算时，我们需要对参数进行调节和优化。
具体来说，我们可以通过多
次求解不同的二阶抛物型偏微分方程，来不断调节参数并逐步优化计
算结果。
以上是一个含参数的高精度两层差分格式，可以用于解决二阶抛物型
偏微分方程的计算问题。
该方法能够保证数值计算的高精度和准确性，同时也能够应对实际问题中的参数不确定性。
李荣华二维抛物方程的初边值问题在数学和工程领域中扮演着重要角色。
通过使用Matlab，我们能够更深入地理解和解决这一问题。
在本文中，我们将从基础概念开始，逐步深入讨论李荣华二维抛物方程的初边值问题，并结合Matlab进行实际分析和解决。
一、初边值问题的概念和涵义1. 初边值问题的定义初边值问题是指在偏微分方程中，除了部分边界条件外，还需给出一些初始条件。
李荣华二维抛物方程的初边值问题即是在方程中给定初值条件和边界条件的问题。
2. 李荣华二维抛物方程简介李荣华二维抛物方程是描述热传导、扩散等现象的数学模型，在物理学和工程领域有着广泛的应用。
它的形式通常为一个关于未知函数和其在空间和时间上的导数的方程。
3. Matlab在初边值问题中的应用Matlab作为一种强大的数学建模和仿真工具，能够帮助我们求解各种类型的偏微分方程，包括初边值问题。
通过Matlab，我们可以对李荣华二维抛物方程进行数值求解和模拟，从而更清晰地理解问题本质。
二、李荣华二维抛物方程的初边值问题解析1. 方程的建立我们需要建立李荣华二维抛物方程的数学模型，并给定相应的初值条件和边界条件。
在Matlab中，我们可以通过定义矩阵和方程的形式来实现这一过程。
2. 数值求解我们利用Matlab提供的数值求解方法，如有限差分法或有限元法，对初边值问题进行求解。
通过程序的编写和运行，我们可以得到方程的数值解，并对物理过程有更直观的认识。
3. 结果分析在求得数值解之后，我们可以对结果进行分析和可视化。
通过Matlab 的绘图功能，我们能够直观地观察李荣华二维抛物方程的演化过程，以及初值条件和边界条件对解的影响。
三、个人观点和总结在本文中，我们深入探讨了李荣华二维抛物方程的初边值问题，并结合Matlab进行了实际分析和解决。
通过对问题的全面评估和数值求解，我们对方程的特性和行为有了更深入的理解。
个人观点上，我认为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 取值。
假定第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.3713； 2-数：0.1447； ∞-数：0.6086（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);%********************************************************************。
