关于解二维椭圆型高精度差分方程的交替方向迭代法

合集下载

二维波动方程的高精度交替方向隐式方法

二维波动方程的高精度交替方向隐式方法

二维波动方程的高精度交替方向隐式方法马月珍;李小纲;葛永斌【摘要】基于二阶微商的四阶紧致差商逼近公式及加权平均思想,提出了数值求解二维波动方程的2种精度分别为O(τ~2+h~4)和O(τ~4+h~4)的交替方向隐式(ADI)格式,以及与其相匹配的第一个时间层的同阶离散格式,并且通过Fourier方法分析了格式的稳定性.该方法在沿每个空间方向上只涉及3个网格基架点,因此可以重复采用TDMA算法,从而大大节省计算时间.数值实验验证了所用方法的精确性和可靠性.【期刊名称】《四川师范大学学报(自然科学版)》【年(卷),期】2010(033)002【总页数】5页(P179-183)【关键词】波动方程;高阶紧致格式;交替方向隐式方法;稳定性【作者】马月珍;李小纲;葛永斌【作者单位】宁夏大学,应用数学与力学研究所,宁夏,银川,750021;宁夏大学,应用数学与力学研究所,宁夏,银川,750021;宁夏大学,应用数学与力学研究所,宁夏,银川,750021【正文语种】中文【中图分类】O241.82有限差分法[1-3]是数值求解偏微分方程的常用方法之一.在航空、气象、海洋、水利等许多流体力学的问题中,常常遇到双曲型偏微分方程,针对传统的差分离散格式,普遍有着精度低且受到很强的稳定性条件限制的缺陷,因此发展其高精度且稳定性好的差分离散方法具有十分重要的意义[4-7].交替方向隐式(AD I)方法是数值求解该类问题的一种非常有效的方法,它将高维问题转化为若干一维问题进行求解,而一维问题又可采用高效的 TDMA算法,从而可以大大提高计算效率,节省存储空间.传统的AD I方法如 Peaceman-Rachford(P-R)格式和Beam-War ming[8]格式都是二阶精度的.另一方面,文献[9]提出了求解二维非定常对流扩散方程的高精度AD I格式,文献[10]提出了求解高维热传导方程的高精度AD I格式.本文依然根据AD I方法的基本思想,提出两种求解二维波动方程的高精度紧致AD I差分格式,为此考虑初边值问题:其中Ω ={(x,y):0≤x,y≤1},Γ为Ω的边界,u(x, y,t)为待求未知量,f(x,y,t)为源项,φ(x,y),ψ(x,y)和 g(x,y,t)均为已知函数且具有充分的光滑性.1 高阶紧致AD I格式用τ表示时间步长,空间取等间距网格,步长用h表示.网格点为(xi,yj,tn),xi=ih,yj=jh,tn = nτ,i,j=0,1,…,N,h=1/N,n≥0.1.1 AD I(2,4)格式考虑(1)式在 n时刻值,对时间导数项采用中心差分,空间导数项采用Kreiss[11]提出的四阶紧致差分公式:则有其中在空间方向以和的算术平均值代替可得对(5)式进行整理且略去高阶项可得为了构造AD I差分格式,采用与文献[9-10]类似的技巧,在(6)式左端加上可得显然,(7)式与(6)式的截断误差同阶,利用可将(7)式写为如下形式引入一个过渡变量则可将 (8)式写为对于过渡变量的边界条件,可以由下式给出(9)式即为求解二维波动方程的高精度紧致 AD I格式,其精度为O(τ2+h4),记为AD I(2,4).1.2 AD I(4,4)格式对时间和空间导数项均采用四阶紧致差分公式,考虑(1)式在 n时刻值,可得对上式进行整理且略去高阶项可得可将(10)式写为如下形式为了构造AD I差分格式,采用与文献[9-10]类似的技巧,在 (11)式左端加上并进行因式分解,可得显然,(11)式与(12)式的截断误差同阶,引入一个过渡变量则 (12)式可写为对于过渡变量的边界条件,可以由下式给出(13)式即为求解二维波动方程的高精度紧致 AD I格式,其精度为O(τ4+h4),记为AD I(4,4).1.3 初始条件的离散因为格式是 3层的.即每一次时间推进都需要知道前两个时间步的值,初始时刻有(2)式精确给出,第一个时间步的值由 (3)式给出,因此,须对 (3)式进行离散,下面推出与(9)和 (13)式相匹配的第一个时间步的离散格式.利用 Taylor展开式将在处展开可得利用(1)~(3)式,且略去高阶项,即可得与(9)式相匹配的第一个时间步的离散格式与(13)式相匹配的第一个时间步的离散格式2 稳定性分析下面,采用 Fourier分析方法对格式进行稳定性分析.引理[12] 实系数二次方程λ2-bλ-c=0的根按模不大于 1的充要条件为|b|≤1-c≤2.2.1 AD I(2,4)格式的稳定性定理 1 格式(9)是无条件稳定的.证明显然,格式(9)可写成用表示采用上述格式进行计算产生的误差,设源项 f无误差存在,则格式的误差项满足格式相应的齐次方程,即记其特征项表示为其中为虚数单位,ηn为第 n个时间层上的波幅,σ1、σ2为波数,令λ=τ/h,则可得误差的传播矩阵为其中矩阵的特征方程为所以有由于Px≥0、Py≥0、Qx≥0、Qy≥0,故|b|≤2,由引理可得格式(9)是无条件稳定的.2.2 AD I(4,4)格式的稳定性定理 2 格式(13)是条件稳定的,其稳定性条件为证明由相同的分析过程可得格式(13)误差的传播矩阵为矩阵的特征方程为由引理得即上式右端显然成立,考察左端可得解之可得|a|λ≤1/3,即格式(13)是条件稳定的,其稳定性条件为|a|λ≤1/3.3 数值验证对于问题(1)~(4),令问题的精确解为计算是用 Fortran77语言进行编程且在 Pentium IV/2.4G PC机上双精度制下进行的.由于2种格式所得线性方程组均为三对角线型,所以对每一步可以采用TDMA 算法.表 1给出了不同网格步长下,问题在 t=0.125时刻,当τ=h时二阶AD I格式[8]、FULL(4,4)格式[7]和τ=h2时AD I(2,4)格式的数值计算结果的最大误差 E和收敛阶 rate=ln(E1/E2)/ln2(E1和E2分别为粗网格及相邻的细网格上的最大误差).结果表明,3种格式均达到了各自的精度,并且AD I(2,4)格式的计算结果要比其它两种格式精确得多.表 2和表 3给出了不同网格步长下,问题在 t =0.25时刻,当τ =h/2时FULL(4,4)格式[7]、AD I(4,4)格式和τ=h2时 HWAL(2,4)[7]格式的数值计算结果的最大误差 E、收敛阶 rate和 CPU时间.结果表明,3种格式均达到了各自的精度,并且AD I(4,4)格式的计算结果要比其它两种格式精确得多,并且计算时间最少.表 4给出了 h=1/64时,问题在 t=3.2时刻,对不同网格比λ(此时τ=λh不同),3种AD I 格式的收敛性的比较.结果表明,当时,AD I(4,4)格式是发散的,即条件稳定的,其它格式仍然是收敛的,即无条件稳定的.这与本文的理论分析结果一致.本文的方法可推广到三维波动方程的数值求解中去,将另文报道.表 1 t=0.125时刻不同空间步长 h的最大误差和收敛阶Table 1 Max imum error and convergencerate att=0.125 for differenthh 二阶ADI格式[8]FULL(4,4)格式[7]ADI(2,4)格式E rate E rate E rate 1/8 2.81e-2 1.29e-34.81e-4 1/16 7.57e-3 1.89 8.20e-5 3.98 3.01e-5 3.97 1/32 1.93e-3 1.975.15e-6 3.99 1.88e-6 4.00 1/64 4.83e-4 1.99 3.22e-7 4.00 1.18e-7 4.001/128 1.22e-4 2.00 1.82e-8 4.14 7.35e-9 4.00表 2 t=0.25时刻不同空间步长 h的最大误差和收敛阶Table 2 Max imum error and convergencerate att=0.25 for differenthh HWAL(2,4)格式[7]FULL(4,4)格式[7]ADI(4,4)格式E rate E rate E rate 1/8 1.54e-2 2.47e-4 2.87e-5 1/162.46e-5 9.29 1.56e-53.98 1.73e-64.05 1/32 1.53e-6 4.00 9.75e-7 4.001.07e-7 4.01 1/64 8.96e-8 4.09 5.99e-8 4.02 6.68e-9 4.00 1/128 5.60e-9 4.002.73e-9 4.45 4.17e-10 4.00表 3 t=0.25时刻不同空间步长 h的 CPU时间Table 3 CPU t ime att=0.25 for differenthλ HWAL(2,4)格式[7]FULL(4,4)格式[7] ADI(4,4)格式1/8 <10e-8 <10e-8 <10e-8 1/16 0.125 0.015 0.015 1/32 1.985 0.125 0.031 1/64 57.50 1.391 0.219 1/128 2830.51 28.39 1.687表 4 当 h=1/64,t=3.2时刻,对不同λ的最大误差Table 4 Max imum erroratt=3.2 andh=1/64 for differentλλ 二阶格式[8] ADI(2,4)格式 ADI(4,4)格式0.4 6.78e-4 4.68e-4 1.93e-8 0.8 2.04e-3 1.84e-3 1.49e+51 1.6 7.03e-3 6.87e-3 3.09e+65 3.2 1.96e-2 1.95e-2 2.80e+39参考文献[1]李胜坤,冯民富,李珊.Benjamin-Bona-Mahony方程的有限差分近似解[J].四川师范大学学报:自然科学版,2003, 26(4):363-365.[2]吕胜关.一类高阶方程组的差分方法[J].郑州大学学报:理学版,1999,32(1):12-19.[3]罗明英,舒国皓,王殿志.RLW方程的有限差分逼近[J].四川师范大学学报:自然科学版,2001,24(2):138-143.[4]Visher J,Wandzura S,White A.Stable,high-order discretization forevolution of the wave equation in 1+1 dimensions[J].J ComputPhys,2004,194:395-408.[5]Wandzura S.Stable,high-order discretization for evolution of the wave equation in 2+1 dimensions[J].J Comput Phys,2004, 199:763-775.[6]葛永斌,吴文权,田振夫.二维波动方程加权平均隐格式及其多重网格方法[J].上海理工大学学报,2002,24(3):205-208.[7]葛永斌,吴文权,田振夫.二维波动方程高精度隐式格式及其多重网格方法 [J].厦门大学学报:自然科学版,2003, 42(6):691-696.[8]孙志忠.偏微分方程数值解法[M].北京:科学出版社,2005.[9]Karaa S,Zhang J.High orderADImethod for solving unsteady convection-diffusion problem[J].J Comput Phys,2004,108:1-9.[10]葛永斌,田振夫,吴文权.高维热传导方程的高精度交替方向隐式方法[J].上海理工大学学报,2007,29(1):55-58.[11]Hirsh R S.Higher order accurate difference solutions of fluid mechanics problem by a compact differencing technique[J].J Comput Phys,1975,19:90-109.[12]陆金甫,关冶.偏微分方程数值解法[M].北京:清华大学出版社,1987.。

二维波动方程的高精度交替方向隐式方法

二维波动方程的高精度交替方向隐式方法

二维波动方程的高精度交替方向隐式方法二维波动方程是一种常见的偏微分方程,用于描述二维空间中的波动现象。

在数值求解二维波动方程时,高精度交替方向隐式方法是一种常用的数值方法。

这种方法在时间和空间方向上交替进行计算,通过隐式格式来增加算法的稳定性和精度。

本文将详细介绍二维波动方程的高精度交替方向隐式方法。

首先,我们来定义二维波动方程的数学模型。

假设二维波动方程的解为u(x,y,t),那么数学模型可以表示为:∂²u/∂t²=c²(∂²u/∂x²+∂²u/∂y²)其中,c表示波速,x,y表示空间坐标,t表示时间。

根据数学模型,我们可以推导出算法的差分格式。

对时间求二阶导数,我们可以使用中心差分公式,得到:(∂²u/∂t²)_(i,j)≈[u(i,j,t+Δt)-2u(i,j,t)+u(i,j,t-Δt)]/Δt²对空间求二阶导数,我们可以使用中心差分公式,得到:(∂²u/∂x²)_(i,j)≈[u(i+1,j,t)-2u(i,j,t)+u(i-1,j,t)]/Δx²(∂²u/∂y²)_(i,j)≈[u(i,j+1,t)-2u(i,j,t)+u(i,j-1,t)]/Δy²将上述差分格式代入数学模型中,得到:[u(i,j,t+Δt)-2u(i,j,t)+u(i,j,t-Δt)]/Δt²=c²{[u(i+1,j,t)-2u(i,j,t)+u(i-1,j,t)]/Δx²+[u(i,j+1,t)-2u(i,j,t)+u(i,j-1,t)]/Δy²}对于高精度交替方向隐式方法,时间和空间的差分方程交替进行迭代,以提高算法的稳定性和精度。

首先,我们可以使用空间差分方程来迭代时间。

时间t的取值范围为[0,T],每个时间步长为Δt。

二维抛物线方程数值解法(ADI隐式交替法)方法

二维抛物线方程数值解法(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。

谱方法和边界值法求解二维薛定谔方程硕士学位论文 精品推荐

谱方法和边界值法求解二维薛定谔方程硕士学位论文 精品推荐

摘要薛定谔方程是物理系统中量子力学的基础方程,它可以清楚地说明量子在系统中随时间变化的规律。

通过求解微观系统所对应的薛定谔方程,我们能够得到其波函数以及对应的能量,从而计算粒子的分布概率,进一步来了解其性质。

在化学和物理等诸多科学研究领域当中,薛定谔方程求解的结果都与实际很相符。

近年来,很多学者通过各种方法研究具有复杂势函数的薛定谔方程,解释了很多重要的物理现象,因此对薛定谔方程的求解具有相当重要的意义。

本文主要是用Galerkin-Chebyshev谱方法和边界值法求解二维薛定谔方程。

首先运用Galerkin-Chebyshev谱方法来对空间导数进行近似,离散二维薛定谔方程,从而将原问题转化为复数域上的线性常微分方程组。

然后用边界值法求解该方程组,所求得的数值解即为原问题的解,之后进行误差分析。

最后利用Matlab进行数值模拟,给出数值解的图像以及误差曲面图像,结果显示此方法精度高且具有很好的稳定性。

关键词:薛定谔方程;Galerkin-Chebyshev谱方法;边界值法;数值解;精度高;稳定AbstractThe Schrödinger equation is the basic equations of quantum mechanics in the physical system. It can clearly describe the regular of the quantum evolves over time. By solving the Schrödinger equation which the micro system correspond, we can get the wave function and energy, and thus calculate the probability distribution of the particles, further understand the nature of it.In chemistry, physics and other fields of scientific research, the results of solving the Schrodinger equation are basically consistent with the actual.In recent years, many researchers used a variety of methods to investigate the Schrödinger equation with complex potential function, and explained a lot of important phenomena.Thus solving the Schrödinger equation has very important significance.The main purpose of this paper is to solve the two dimensional Schrödinger equation through the Galerkin-Chebyshev spectral method and the boundary value method. First we use the spectral method to approximate the spatial derivation, discretize the two dimensional Schrödinger equation,and transform the original problem into a set of linear ordinary differential equations in the complex number field.Then by using the boundary value method to solve the equations, that the numerical solutions is the solutions of the original problem, and then analyze the error. Finally we use Matlab to conduct the numerical simulation, and give the images of the numerical solutions and errors, which show that the methods have high precision and good stability.Keywords: Schrödinger equation, Galerkin-Chebyshev spectral method, boundary value method, numerical solutions, high precision, stability目录摘要 (I)Abstract (II)第1章绪论 (4)1.1课题研究的背景和意义 (4)1.2国内外研究现状 (5)1.3本文的主要研究内容 (5)第2章预备知识 (7)2.1克罗内克积的简介 (7)2.2Chebyshev多项式介绍及其性质 (8)2.3Chebyshev正交逼近的性质 (9)2.4投影算子的性质 (10)2.5本章小结 (11)第3章Galerkin-Chebyshev谱方法和边界值法 (12)3.1用Galerkin-Chebyshev谱方法求解椭圆型方程 (12)3.2用边界值法求解常微分方程 (13)3.3本章小结 (17)第4章求解二维薛定谔方程 (18)4.1区域和边界条件的处理 (18)4.1.1 区域的处理 (18)4.1.2 边界条件的处理 (20)4.2二维薛定谔方程的求解 (23)4.3误差分析 (24)4.4本章小结 (29)第5章数值模拟 (30)结论 (35)参考文献 (36)哈尔滨工业大学学位论文原创性声明及使用授权说明 .....错误!未定义书签。

二维hammerstein方程数值解的高精度算法

二维hammerstein方程数值解的高精度算法

二维hammerstein方程数值解的高精度算法二维Hammerstein方程是一种常见的非线性微分方程,由线性动力系统和非线性静态非连续增量函数组成。

它在信号处理、图像处理、控制理论等领域具有广泛的应用。

解决二维Hammerstein方程的数值方法一直是学术界的热点问题之一,因为它具有挑战性和重要性。

在本篇文章中,我们将讨论关于二维Hammerstein方程数值解的高精度算法。

1. 了解二维Hammerstein方程我们需要了解二维Hammerstein方程的基本形式和特点。

该方程可以被表示为:\[ y(n, m) = \int_{0}^{1} \int_{0}^{1} K(n, m, x, y)f(x, y)dxdy \]其中,\(y(n, m)\)代表方程的输出,\(f(x, y)\)是方程的输入,而\(K(n, m, x, y)\)是非线性增量函数。

二维Hammerstein方程描述了输入和输出之间的关系,通过非线性增量函数将输入信号映射到输出信号。

2. 数值解的常见方法对于二维Hammerstein方程的数值解,目前存在许多常见方法。

其中最常见的方法是使用离散化的方式来近似方程中的积分。

通过将定义域离散化为多个细小的网格,并使用数值积分方法(如Simpson法则或梯形法则),可以得到方程的数值解。

这些方法相对简单,但在某些情况下可能会导致精度降低。

3. 高精度算法为了提高二维Hammerstein方程数值解的精度,我们可以应用一些高精度算法。

这些算法可以通过减小离散化步长、使用更准确的数值积分方法或采用更高阶的数值差分格式来实现。

一种常见的高精度算法是通过使用多项式插值来近似非线性增量函数。

通过在每个离散点处计算非线性增量函数的值,并使用插值方法来估计其他任意点的函数值,我们可以得到更准确的数值解。

其中,拉格朗日插值和样条插值是常用的插值方法,在二维Hammerstein方程的求解中均有应用。

差分方程及其Z变换法求解

差分方程及其Z变换法求解

由图:x1 (k 1)T x2 (kT )
zX 1 ( z ) zx1 (0) X 2 ( z )
x2(kT)
z
1
x1(kT)
z 1
x1(0) 1
x1 ( z)
x2(z) y[(k+1)T]
例2:画出例2所示离散系统的模拟图
y[(k 1)T ] -( KT -1) y(kT ) + KTr(k 1)T y (kT ) T
(T 很小)
(2)
式中:T为采样周期,(2)代入(1)得:
y (k 1)T (KT 1) y(kT ) KTr(kT )
y(k 1) ( K 1) y(k ) Kr (k )
(3)
二、离散系统差分方程的模拟图
连续系统采用积分器s-1作为模拟连续系统微分方程的主要器件; 与此相对应,在离散系统中,采用单位延迟器z-1。 单位延迟器:把输入信号延迟一个采样周期T秒或延迟1拍。
y (t ) y(kT )(t kT )
* n 0
[1/ 6 1/ 2(1) k 2 / 3(2) k ](t kT )
n 0

例4:求y[(k+2)T]+y[(k+1)]+0.24y(kT)=u(kT)在单位阶跃函数 作用下的解。初始条件y(0)=0, y(T)=1。
U ( z ) Z [1(t )]
z z 1
0.446 1.429 1.875 z ( z -1) ( z 0.4) ( z 0.6)
y(kT ) 0.446 1.429(-0.4)k -1.875(-0.6)k
y (t ) [0.446 1.429(0.4) k 1.875(0.6) k ](t kT )

二阶椭圆偏微分方程实例求解(附matlab代码)(完整资料).doc

二阶椭圆偏微分方程实例求解(附matlab代码)(完整资料).doc

【最新整理,下载后即可编辑】《微分方程数值解法》期中作业实验报告二阶椭圆偏微分方程第一边值问题姓名:学号:班级:2013年11月19日二阶椭圆偏微分方程第一边值问题摘要对于解二阶椭圆偏微分方程第一边值问题,课本上已经给出了相应的差分方程。

而留给我的难题就是把差分方程组表示成系数矩阵的形式,以及对系数进行赋值。

解决完这个问题之后,我在利用matlab 解线性方程组时,又出现“out of memory ”的问题。

因为99*99阶的矩阵太大,超出了分配给matlab 的使用内存。

退而求其次,当n=10,h=1/10或n=70,h=1/70时,我都得出了很好的计算结果。

然而在解线性方程组时,无论是LU 分解法或高斯消去法,还是gauseidel 迭代法,都能达到很高的精度。

关键字:二阶椭圆偏微分方程 差分方程 out of memory LU 分解 高斯消去法 gauseidel 迭代法一、题目重述解微分方程:()()2222((,))((,))()(,)()(,)(,)1y x x x y y x y y x xy xy e u x y e u x y x y u x y x y u x y u x y y e x e e y x e --+++-+=-++++已知边界:(0,)1,(1,),(,0)1,(,1)y x u y u y e u x u x e求数值解, 把区域[0,1][0,1]G 分成121/100,1/100h h ,n =100 注:老师你给的题F 好像写错了,应该把22x y y e x e +改成22y x y e x e +。

二、问题分析与模型建立2.1微分方程上的符号说明A (A ,A )=A A A (A ,A )=A A A (A ,A )=A +A A (A ,A )=A −AA (A ,A )=1 A (A ,A )=()()22221y x xy xy y e x e e y x e -++++2.2课本上差分方程的缺陷课本上的差分方程为:A AA A AA−(A A−1,A A A−1,A+A A,A−1A A,A−1+A A+1,A A A+1,A+A A,A+1A A,A+1)=A AA{A A−1,A=A−2(A A−1/2,A+AA AA2)A A,A−1=A−2(A A,A−1/2+AA AA2)A A+1,A=A−2(A A+1/2,A−AA AA2)A A,A+1=A−2(A A,A+1/2−AA AA2)A AA=A−2(A A+1/2,A+A A−1/2,A+A A,A−1/2+A A,A+1/2)+A AA 举一个例子:当i=2,j=3时,A AA=A23;当i=3,j=3时,AA−1,A=A23。

二阶椭圆偏微分方程实例求解(附matlab代码)

二阶椭圆偏微分方程实例求解(附matlab代码)

《微分方程数值解法》期中作业实验报告二阶椭圆偏微分方程第一边值问题姓名:学号:班级:2013年11月19日二阶椭圆偏微分方程第一边值问题摘要对于解二阶椭圆偏微分方程第一边值问题.课本上已经给出了相应的差分方程。

而留给我的难题就是把差分方程组表示成系数矩阵的形式.以及对系数进行赋值。

解决完这个问题之后.我在利用matlab 解线性方程组时.又出现“out of memory ”的问题。

因为99*99阶的矩阵太大.超出了分配给matlab 的使用内存。

退而求其次.当n=10.h=1/10或n=70.h=1/70时.我都得出了很好的计算结果。

然而在解线性方程组时.无论是LU 分解法或高斯消去法.还是gauseidel 迭代法.都能达到很高的精度。

关键字:二阶椭圆偏微分方程 差分方程 out of memory LU 分解 高斯消去法 gauseidel 迭代法一、题目重述解微分方程:()()2222((,))((,))()(,)()(,)(,)1y x x x y y x y yxxyxye u x y e u x y x y u x y x y u x y u x y y e x e e y x e--+++-+=-++++已知边界:(0,)1,(1,),(,0)1,(,1)y x u y u y e u x u x e ====求数值解, 把区域[0,1][0,1]G =?分成121/100,1/100h h ==.n =100 注:老师你给的题F 好像写错了.应该把22x y y e x e +改成22y x y e x e +。

二、问题分析与模型建立2.1微分方程上的符号说明()()22221y x xy xy y e x e e y x e -++++2.2课本上差分方程的缺陷课本上的差分方程为:举一个例子:当i=2,j=3时.;当i=3,j=3时.。

但是.显然这两个不是同一个数.其大小也不相等。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

关于解二维椭圆型高精度差分方程的交替方向迭代法交替方向迭代法是一种常用于解二维椭圆型差分方程的方法。

在许多实际应用中,需要对这种方程进行高精度的求解。

交替方向迭代法可以有效地解决这个问题,并且在计算效率和内存利用方面都具有优势。

交替方向迭代法的基本思想是将二维椭圆型差分方程化为两个一维的差分方程,然后分别对其进行迭代求解。

具体地,将二维椭圆型差分方程表示为:
$-\frac{\partial^2 u}{\partial 某^2}-\frac{\partial^2
u}{\partial y^2}=f(某,y)$。

其中,$u(某,y)$是要求解的未知函数,$f(某,y)$是已知的函数。

接下来,我们可以将上式分别用某轴和y轴方向进行差分,得到两个一维的差分方程:
$-\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{\Delta 某^2}-
\frac{u_{i,j-1}-2u_{i,j}+u_{i,j+1}}{\Delta y^2}=f_{i,j}$。

$-\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{\Delta 某^2}-
\frac{u_{i,j-1}-2u_{i,j}+u_{i,j+1}}{\Delta y^2}=f_{i,j}$。

其中,$u_{i,j}$表示在坐标$(i,j)$处的未知函数值,$\Delta 某$和$\Delta y$分别表示在某轴和y轴方向的差分间隔。

接下来就是交替方向迭代法的核心部分。

首先,我们需要对第一个方程中的某方向进行求解,可以通过迭代求解来不断逼近最终的解。

假设我们用$u^{k}_{i,j}$表示在第k次迭代时,在坐标$(i,j)$处的近似解。

那么对于每个坐标$(i,j)$,我们可以通过以下公式来更新近似解:
$u^{k+1}_{i,j}=\frac{1}{2}(\frac{u^{k}_{i-
1,j}+u^{k}_{i+1,j}}{\Delta 某^2}+\frac{u^{k}_{i,j-
1}+u^{k+1}_{i,j+1}}{\Delta y^2}-f_{i,j})$。

同样的,我们可以用类似的公式对y方向进行迭代求解,直到达到一定的精度或者迭代次数到达预设的阈值为止。

值得注意的是,交替方向迭代法在计算效率和内存利用方面都具有优势,因为它可以将原本的二维椭圆型差分方程化为两个一维的差分方程,并且只需要存储每个坐标处的近似解即可。

因此,可以对大规模的二维椭圆型差分方程进行高效的求解。

总之,交替方向迭代法是一种有效的解二维椭圆型差分方程的方法。

它能够对高精度的差分方程进行求解,并且在计算效率和内存利用方面都具有优势。

相关文档
最新文档