关于解二维椭圆型高精度差分方程的交替方向迭代法
二维波动方程的高精度交替方向隐式方法
二维波动方程的高精度交替方向隐式方法马月珍;李小纲;葛永斌【摘要】基于二阶微商的四阶紧致差商逼近公式及加权平均思想,提出了数值求解二维波动方程的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隐式交替法三种解法及误差分析(一般的教材上只说第一种)理论部分参看孙志忠:偏微分方程数值解法注意: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方程数值解的高精度算法。
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变换法求解
由图: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
【最新整理,下载后即可编辑】《微分方程数值解法》期中作业实验报告二阶椭圆偏微分方程第一边值问题姓名:学号:班级: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代码)
《微分方程数值解法》期中作业实验报告二阶椭圆偏微分方程第一边值问题姓名:学号:班级: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时.。
但是.显然这两个不是同一个数.其大小也不相等。
关于解二维椭圆型高精度差分方程的交替方向迭代法
关于解二维椭圆型高精度差分方程的交替方向迭代法交替方向迭代法(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方法是一种隐式迭代方法,对稳定性和收敛性有较好的保证。
二维波动方程的高精度交替方向隐式方法
二维波动方程的高精度交替方向隐式方法
1二维波动方程
二维波动方程是水文学研究与应用领域中十分重要的数学模型,在流体运动和水利、水文工程中都有广泛的应用。
它解释了水在介质中的波动行为,包含了二维水声传播模型和正弦波运动的数学表达式。
2高精度交替方向隐式方法
上述二维波动方程的微分方程本身被认为是难以求解的非线性方程,而高精度交替方向隐式方法则是这种难以求解的模型的关键解决方案。
简略的说,这种方式的求解针对一般形式的波动方程,分别以x 和y方向上的精度去求解,直到到达满意的解析结果。
3此方法的优点
这种高精度交替方向隐式方法有助于改善现有二维波动方程数值求解结果的准确性,并为在多维情况下求解非线性微分方程提供了坚实的基础。
此外,这种方法可以有效节省更多的计算成本,减少计算时间,为水文学研究与应用领域带来众多的便利。
4结论
总而言之,高精度交替方向隐式方法在求解二维波动方程方面实现了非线性微分方程的准确分析,且有效减少计算成本,其计算结果
也更为精确,为水文学研究与应用领域的应用提供了更便捷的计算方法。
求解非线性方程的三种新的迭代法
求解非线性方程的三种新的迭代法求解非线性方程是数值分析中的一个重要问题,对于复杂的方程往往无法通过解析方法求得精确解。
迭代法成为了求解非线性方程的常用方法之一。
传统的迭代法如牛顿法、割线法等在实际应用中存在着一些问题,近年来研究者提出了一些新的迭代法,以解决传统迭代法的一些缺点。
下面将介绍三种新的迭代法:拟牛顿法、混合法和基因迭代法。
拟牛顿法是一种通过构建迭代更新的矩阵逼近方向导数的方法。
传统的牛顿法每次需要求解并存储Hessian矩阵的逆矩阵,计算复杂度较高。
而拟牛顿法通过不直接计算Hessian矩阵的逆矩阵,而是通过一些迭代更新矩阵来逼近Hessian矩阵的逆矩阵,从而减少了计算复杂度。
其基本思想是通过更新矩阵不断逼近Hessian矩阵的逆矩阵,使得迭代得到的解更加接近最优解。
拟牛顿法在求解非线性方程时具有收敛速度快、计算精度高等优点,在实际应用中得到了广泛的应用。
另一种新的迭代法是混合法,混合法是将不同的迭代方法进行组合,通过不同的方式交替使用各种迭代方法,以期望达到更好的收敛效果。
将牛顿法和割线法进行混合,每次迭代交替使用两种方法,通过对两种迭代方法的优势互相弥补,从而提高收敛速度和稳定性。
混合法在实际应用中表现出了良好的性能,特别是在求解高阶非线性方程或者函数具有复杂结构的情况下,混合法能够更快速地收敛到最优解。
最后一种新的迭代法是基因迭代法,基因迭代法是受到生物进化理论的启发而提出的一种迭代方法。
其基本思想是通过模拟生物进化中的遗传、变异、选择等过程,不断优化迭代的解。
在每一代迭代中,通过交叉、变异等操作产生新的解,并通过适应度函数进行选择,使得优秀的解得以保留并得到进化,从而逐步优化迭代的解。
基因迭代法在求解非线性方程时可以克服局部最优解的问题,找到更接近全局最优解的解。
在求解复杂的非线性方程时,基因迭代法具有明显的优势。
以上介绍了三种新的迭代法:拟牛顿法、混合法和基因迭代法。
这些新的迭代法在求解非线性方程时有着良好的性能,能够更快速地收敛到最优解,具有更好的稳定性和精度。
差分方程的迭代解法
1
实验目的
学习用Matlab计算离散信号的功率和能量。
学习并掌握用迭代法求解差分方程的方法 。 掌握用Matlab进行离散卷积运算的数值方法 和解析方法。加深对离散卷积的理解。
2
实验原理与说明
离散信号的能量与功率
与连续信号类似,离散信号也可分为能量信号和 功率信号。对于非周期信号,信号能量定义为
1
0
其它
2
0
其它
解 用Matlab并调用DSCONV()函数,程序如下:
% 计算离散信号的卷积 exp23_2a.m n1=-2;f1=[2 2 2 2 2 2]; % 序列的起始点,序列值 n2=0;f2=[1 1 1 1 1 1]; % 序列的起始点,序列值 M=6; % 将卷积值显示在中间,左右插入M点 dsconv(f1,n1,f2,n2,M) 在命令窗口显示的卷积结果 y= 2 4 6 8 10 12 10 8 6 4 2
9
计算示例
例2 求下述差分方程的解 y(k ) 1.5 y(k 1) y(k 2) 2 f (k 2) y (2) 2 。 其中输入信号 f (k ) (k ) ,初始条件 y (1) 1 ,
解 Matlab程序如下: % 计算例2的程序 exp23_1.m a=[-1.5 1];b=[0 0 2]; y0=[2 1];f0=[0 0]; n=0:30; f=ones(1,length(k)); y=recur(a,b,n,f,f0,y0); stem(n,y,'.'), xlabel('k'),ylabel('y(k)')
第二项与上类似
椭圆方程迭代法介绍
第三章 椭圆型问题的差分法§3-1 流体力学中的椭圆型问题·无旋流场中 速度势02=∇ϕ (Laplace Eq.) ·二维不可压定常流动,利用涡-流函数表示:⎪⎩⎪⎨⎧-=∇∇=∇⋅+∂∂方程Poisson V tζϕζυζζ22·不可压分离流问题中,扰动压力场:02='∇p·定常的N -S 方程求解问题·在网格自动生成中,求解椭圆型方程的网格生成方法 由于椭圆型方程的数学性质:求解域内部任何一点的解函数依赖于所有边界上的边界条件,因此从数值计算方法来看,就不能从一部分边界起步进行推进计算到另外的边界,这与发展方程的求解方法有很大的差别,椭圆型方程的数值求解方法,只能是在整个流场中进行迭代计算来求解。
§3-2 椭圆型问题的迭代法求解(一)迭代法的基本概念例:方程ζψ=∇2 ( Poisson 方程) 二维ζψψ=∂∂+∂∂2222yx差分离散j i j i j i j i ji j i j i y x ,21,1,1,2,1,,1)(2)(2ζψψψψψψ=∆+-+∆+---+-+ ………………..(*)写成矩阵形式代数方程组为: B A =ψ (1)其中 ⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=A ⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡=ψJ I j i ,,2,11,1ψψψψ 一般地,对于线性方程组有B A =ψ,欲求未知函数ψ的解矢量 若A 为非奇异矩阵,即:∃1-A ,则B A 1-=ψ由于A 是个阶数甚大的矩阵(非三对角),直接求解,或利用Gauss 消去法求逆矩阵,计算量及所需计算机的内存都将十分巨大,所以在实际计算中不希望采用直接法求解。
迭代法的基本思想是:定义一个序列)(k ψ,当∞→k 时,B A k 1)(-→ψ,从而得到方程(1)的解。
迭代法设法给出[])()2()1()(,,,,r k k k k k B A F ---ψψψ=ψ 的迭代关系。
偏微分方程(椭圆型)数值解2.4-5
ai−1, j
=
⎛ ⎜⎝
Ai
−
1 2
,
j
−
h1 2
Cij
⎞ ⎟⎠
h1−2
,
ai, j+1
=
⎛ ⎜⎝
Ai
∂n q1q2
∂n q2q3
∂n q3m4
• = ∫∫ f ( x, y)dxdy m1
G0
p1
(4.6)
对于后四项仿照公式(4.3)的方法离散化,例如
∫ ∫ ( ) m1q1
∂u ds ≈ ∂n
m1q1 ⋅ p0 p1
u1 − u0
,
( ) q2q1
∂u ds ≈ ∂n
q1q2 ⋅ p0 p2
u2 − u0
考虑矩形网格:h1和h2分别为x和y方向的步长,Gh为网格 节点集合,Gh为网格界点集合,Gh = Gh ∪ Γh 。 总假设Gh是联通
的,即 ∀P, P ∈ Gh 必有一串 Pi ∈Gh (i = 1, 2, L, m −1) 可将 P, P
排列成顺序序列:P, P1, P2 , L, Pm , P 使前后两点互为相邻节点。
G0的面积:
m
q7
(G
=
0)
q1 ,i
⎛ = 12× ⎜⎜⎝
= 1,
1× 2
2,
3 3
L, 6。
h
×
1 2
h
⎞ ⎟⎟⎠
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 取值。
椭圆 数值解
椭圆型方程的差分方法[五点格式] 考虑拉普拉斯方程的第一边值问题式中为定义在D的边界S上的已知函数.采用正方形网格,记,在节点(i, j)上分别用差商代替,对应的差分方程为(1)或即任一节点(i, j)上的值等于周围相邻节点上解的值的算术平均,这种形式的差分方程称为五点格式,在边界节点上取(2)式中是与节点(i, j)最接近的S上的点.于是得到了以所有内节点上的值为未知量的若干个线性代数方程,由于每一个节点都可列出一个方程,所以未知量的个数与方程的个数都等于节点的总数,于是,可用通常的方法(如高斯消去法)解此线性代数方程组,但当步长不很大时,用高斯消去法将会遇到很大困难,可用下面介绍的其他方法求解.若h 0时,差分方程的解收敛于微分方程的解,则称差分方程为收敛的. 在计算过程中,由于进行四则运算引起舍入误差,每一步计算的舍入误差都会影响以后的计算结果,如果这种影响所产生的计算偏差可以控制,而不至于随着计算次数的增加而无限增大,则称差分方程是稳定的.[迭代法解差分方程] 在五点格式的差分方程中,任意取一组初值{},只要求它们在边界节点(i , j )上取以已知值,然后用逐次逼近法(也称迭代法)解五点格式:逐次求出{}.当(i+1, j ),(i -1, j ),(i , j -1),(i , j+1)中有一点是边界节点时,每次迭代时,都要在这一点上取最接近的边界点的值.当n →∞时,收敛于差分方程的解,因此n 充分大时,{}可作差分方程的近似解,迭代次数越多,近似解越接近差分方程的解.[用调节余数法求节点上解的近似值] 以差商代替Δu 时,用节点(i+1, j ),(i-1, j ),(i , j+1),(i , j -1)上u 的近似值来表示u 在节点(i , j )的值将产生的误差,称此误差为余数,即设在(i , j )上给以改变量,从上式可见将减少4,而其余含有的差分方程中的余数将增加,多次调整的值就可将余数调整到许可的有效数字的范围内,这样可获得各节点上u (x , y )的近似值.这种方法比较简单,特别在对称区域中计算更简捷.例 求Δu =0在内节点A ,B ,C ,D 上解的近似值.设在边界节点1,2,3,4上分别取值为1,2,3,4(图14.8) 解 记u (A)=,点A ,B ,C ,D 的余数分别为图14.8-4+ u B+ u c +5=-4 u B + u D+7=R B-4 u c+ u D+3=R Cu B+ u c-4u D+5=R D以边界节点的边值的算术平均值作为初次近似值,即=u B(0)=u C(0)=u D(0)=2.5则相应的余数为:=0, R B=2, R C= -2, R D=0最大余数为±2.先用δu C=-0.5把R C缩减为零,u C相应地变为2,这时,R D也同时缩减(-0.5),新余数是=-0.5,R B=2,, R D=-0.5.类似地再变更δu B=0.5,从而u B变为3,则得新余数为.这样便可消去各节点的余数,于是u在各节点的近似值为:=2.5, u B=3, u C=2, u D=2.5现将各次近似值及余数列表如下:[解重调和方程的差分方法] 在矩形D(x0≤x≤x0+a,y0≤y≤y0+a)中考虑重调和方程取步长,引直线族(i, j = 0, 1, 2n)作成一个正方形网格.用差商代替偏导数上式表明了以(x, y)为中心时,u(x, y)的函数值与周围各点函数值的关系,但对于邻近边界节点的点(x, y),如图14.9中的A,就不能直接使用上式,此时将划分网格的直线族延伸,在延伸线上定出与边界距离为h的点,称这些点为外邻边界节点,如图14.9以A为中心时,点E,C为边界节点,点J,K为E,C的外邻边界节点,用下法补充定义外邻边界节点J处函数的近似值u J,便可应用上面的公式.1 边界条件为图14.9时,定义.2︒边界条件为时,定义.[其他与Δu有关的网格]1︒三角网格(图14.10(a))取P0(x, y)为中心,它的周围6个邻近节点分别为:则式中,R表示余项.2︒六角网格(图14.10(b))取P0(x, y)为中心,它的三个邻近节点分别为则.图14.103 极坐标系中的网格(图14.10(c))取P 0(r,)为中心,它的四个邻近节点分别为而拉普拉斯方程的相应的差分方程为。
Newton迭代法求解高阶微积分方程的误差分析
Newton迭代法求解高阶微积分方程的误差分析在实际工程应用领域,存在很多复杂的微积分方程,这类微积分方程难以使用初等的解法进行求解,其解析表达式的找出也十分的困难。
为此,一般使用數值求解法解决这类复杂的微积分方程。
随着电脑技术的发展,高阶微积分方程数值求解方式也得到了一定的发展,该种求解方式能够很好的解决工程技术领域和自然科学领域的复杂方程。
就现阶段来看,高阶微积分数值求解方式已经成为数学工作人员研究的重点问题之一。
下面就针对一类四阶微积分方程的差分迭代解法进行深入的分析。
一、四阶微积分方程模型的建立一桁车横梁模型,其中桁车横梁长度为L,在竖直方向的变形位移与载荷可以用下述方程描述:(1.1)其中;载荷为桁车横梁自重与承载之和。
由桁车的支撑边界,其该方程的初始条件为:。
在式(1.1)中,均为常数,为上的连续非负(或非正)函数,其正负取决于的定义方向。
目前,对此类高阶方程已进行了深入的分析和研究。
本文结合有限差分法,利用Newton迭代的求解方法求解该方程。
下面将该方程求解转化为对非线性方程的求解:式(1.5)为一非线性方程。
这样就将初始微积分模型的求解转化为求解方程式(1.5)。
式(1.5)解的性质取决于函数的性质。
二、Newton迭代法求解方程步骤一般情况下,对于的解析解的求解很复杂。
Newton迭代法在处理非线性方面有着一定的优势,使用牛顿迭代法能够将简化为:再将代入式(2.1)中,求出。
并将和差的绝对值与期望所求精度进行比较,为期望方程所求精度的一个无限小。
如果,那么迭代终止,否则将用代替,转入下一轮继续循环计算。
通俗的说,Newton迭代法其实也就是采用切线逼近的办法求得一个在我们求解公差范围内的一个合理解。
求解精度越小,计算结果的精度越高,但迭代次数会增加。
对的两种不同形式的差商,分别定义为Newton法1和Newton 法2。
三、Newton迭代法计算高阶方程结果分析为验证该高阶方程的迭代求解结果并分析其误差,对式(1.1)中的参数取值为:,;令该微积分方程的解为:,则可计算出,再设定求解期望精度=10-8。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 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方向进行迭代求解,直到达到一定的精度或者迭代次数到达预设的阈值为止。
值得注意的是,交替方向迭代法在计算效率和内存利用方面都具有优势,因为它可以将原本的二维椭圆型差分方程化为两个一维的差分方程,并且只需要存储每个坐标处的近似解即可。
因此,可以对大规模的二维椭圆型差分方程进行高效的求解。
总之,交替方向迭代法是一种有效的解二维椭圆型差分方程的方法。
它能够对高精度的差分方程进行求解,并且在计算效率和内存利用方面都具有优势。