一类二维抛物型方程的ADI格式

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)τ取值。

抛物方程的ADI方法的几点说明摘要:在此文中,我们介绍了抛物方程的几种新的ADI方法。

从这些方法的构建可以看出, ADI法是一个开放式的算法,具有很强的可重构性。

进而启发我们将日常学习到的理论知识举一反三地应用到实际生产中去。

关键词:ADI方法高精度变系数高维方程产生于上世纪五十年代的交替方向隐式(ADI)方法[1,2,3],在数值数学研究与应用领域占有十分重要的地位。

该方法能把多维问题化为一系列的一维问题处理,而且一维问题所要求解的方程组为三对角的,很容易使用熟知的追赶法在计算机上实现。

该方法具有运算速度快、存储量小、无条件稳定等优点,几十年来一直受到人们的青睐。

下面我们将根据抛物方程的不同形式和实际问题的需求,介绍几种新型的交替方向隐格式。

1 传统的ADI法第一个ADI法是Peaceman和Rachford[1]于1955年提出的。

为了显示其方法,考虑如下的二维热传导方程:4 结语在这篇文章中,依据高等教材中交替方向隐格式的基本算法理论,我们介绍了抛物方程几种新型交替方向隐格式。

从这些差分格式的构建不难看出,交替方向隐格式法不但具有运算速度快、存储量小、无条件稳定等优点,而且还是一种开放式的数值算法,具有很强的可重构性,进而启发我们将日常学习到的理论知识举一反三地应用到实际生产中去。

参考文献[1] Peaceman D W,Rachford Jr H H,The numerical solution of parabolic and ellipfic differential equations, J SIAM,3,1955:28-41.[2] 李荣化,刘播,微分方程数值解法[M].4版.北京:高等教育出版社,2009.[3] 胡健伟,汤怀民.微分方程数值方法[M].北京:科学出版社,2000.[4] 闵涛,张海燕,周宏宇,等.二维变系数热传导方程初边值问题的交替方向隐格式[J].西安工业大学学报,2007,27(2).[5] 葛永斌,田振夫,吴文权.高维热传导方程的高精度交替方向隐式方法[J].上海理工大学学报,2007(29):55-58.[6] 朱琳琳.抛物方程的交替方向隐格式[D].河南师范大学,2007.。

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 取值。

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方向）的方程组中求解。

通过交替进行，可以逼近二维抛物线方程的解。

一类二阶二维常系数抛物型方程的参数估计方法抛物型方程是一类很常见的偏微分方程，在数学和工程领域都有着广泛的应用。

通常情况下，我们可以通过已知的数据和条件来求解这类方程的解析解，但是有时候我们并不能得到方程的解析解，这时我们就需要考虑用参数估计的方法来求解方程的解。

在本文中，我们将讨论一类二阶二维常系数抛物型方程的参数估计方法。

这类方程的一般形式可以写为：$$\frac{∂^{2} u}{∂t^{2}} = a\frac{∂^{2} u}{∂x^{2}} +b\frac{∂u}{∂x} + cu$$其中$a,b,c$是常数，$u(x,t)$是变量。

接下来我们将介绍一种基于最小二乘法的参数估计方法，该方法可以帮助我们估计方程中的未知参数。

首先，我们假设我们有一些已知的数据点$(x_{i},t_{i},u_{i}),i=1,2,...,n$，我们可以用数据点构建一个误差函数，使其最小化。

误差函数的一种形式可以写为：$$E(a,b,c) = \sum_{i=1}^{n} (u_{i} - u(x_{i},t_{i}))^{2}$$其中$u(x_{i},t_{i})$是抛物型方程的解析解在点$(x_{i},t_{i})$的取值。

我们的目标是找到$a,b,c$的最优估计值，使得误差函数最小化。

接下来，我们可以用最小二乘法来估计$a,b,c$的值。

最小二乘法的基本思想是通过最小化误差函数来估计参数值。

我们可以将误差函数$E(a,b,c)$对$a,b,c$分别求偏导，并让导数为0，得到方程组：$$\begin{cases}\frac{∂E}{∂a} = 0 \\\frac{∂E}{∂b} = 0 \\\frac{∂E}{∂c} = 0 \\\end{cases}$$解这个方程组，我们就可以得到$a,b,c$的最优估计值。

最后，我们可以用已知的数据点和估计值$a,b,c$来重新求解抛物型方程，得到近似解，并且可以用残差分析方法来检验我们的估计结果是否准确可靠。