隐式有限差分编程

交替方向隐式时域有限差分法中的Berenger理想匹配层

中围分章鳙号：０１２０（０６０ —４８０１０ —４０２０）３０５ —４
Ｂｅｅｅ ’ ｅｆｃｌｔｈｅａｅｏｈｌｅｎｔｎｒｎｇｒｓｐｒｅｔｙｍａｃｄｌｙｒｆｒｔｅａｔｒａｉｇ
ｉｌｉ（ｍｐｉｔＡＤＤｌｏｉｍ；ｅｆｃｌｔｈｄｌｙｒＰＭＬ）ｃａｇｒｔｈｐｒｅｔｙｍａｃｅａｅ（
当利用时域有限差分法（ＤＴ模拟开区域的电磁场问题时，收边界条件（Ｃ）必须的．ｅｅｇｒＦＤ）吸ＡＢ是Ｂｒｎｅ提出的理想匹配层（ＭＬ是一种建立在场分裂基础上的媒质吸收边界条件［，种ＡＢＰ）１这］Ｃ理论上可实现反射误差任意小．无条件稳定的交替方向隐式时域有限差分法（ＩＤＡＤ— ＴＤ）法［出现以后，何实现ＡＤＩＤＦ算２如 — ＴＤ和Ｆ
常的时域有限差分法推广到交替方向隐式时域有限差分法．已有的离散方案比较。和采用文中提出的商
散方案可使理想匹配层吸收边界的反射误差成数量级地减小．关键词：交瞽方向隐式方法ｌ时域有限差分法Ｉ想匹配层理
ｉｓｅｄｏｈｔｎａｄＦＤｎｔａｆｔｅｓａｄｒｌＴＤｔｏ．Ｃｍｐｒｄｗｉｈｉｃｅｉｔｎｓｈｍｅｐｏｏｅｍｅｈｄｏａｅｔｔｅｄｓｒｔａｉｃｅｒｐｓｄｈｚｏｅｒｅ，ｔｅｎｗｎｋｓｔｅｒｆｃｉｎｅｒｒｆＰＬａｓｒｉｇｂｕｄｒｏｄｔｎｄｃｅｓａｌｒｈｅｏｅｍｅｈｅｅｔｒｏｓｏＭｂｏｂｎｏｎａｙｃｎｉｏｅｒｅｉａｌｏｉａ

第五讲——显式差分和隐式差分(5)(格式整齐)

左端：n+1时刻的值；右端：n时刻的值。
特点：结构简洁，直接求解，求解速度快。
但是，时间步长需满足：
显式差分格式才能得到稳定的数值解，否则，数值解将会不稳定而振荡。
高级材料
19
显示差分格式示意图
高级材料
20
2. 隐式差分格式：
高级材料
时间一阶精度空间二阶精度
21
隐式有限差分格式
高级材料
22
初始条件：
高级材料
26
内部节点：
A = sparse(nx,nx); for i=2:nx-1 A(i,i-1) = -s; A(i,i ) = (1+2*s); A(i,i+1) = -s; end
边界节点：
A(1 ,1 ) = 1; A(nx,nx) = 1;
载荷项：
rhs = zeros(nx,1); rhs(2:nx-1) = Told(2:nx-1); rhs(1) = Tleft; rhs(nx) = Tright; 高级材料
end
end
高级材料
b=a^(-1); c=zeros(135,1); for i=121:135
c(i,1)=25;end d=b*c; s=zeros(11,17); for i=2:16
s(11,i)=100; end for i=1:9 for j=1:15; s(i+1,j+1)=d(15*(i-1)+j,1); end end
内部
边界
27
Crank-Nicolson 隐式差分格式的程序实现

sTi
n1 1

(2

2s)Ti n 1

有限差分法

有限差分法一、单变量函数：用中心差分法（matlab程序见附录）计算结果如下：图1 中心差分法表1 数据对比二、一维热传导：在此取φ(x)=0，g1(t)= g2(t)=100-100*exp(-t)问题描述：已知厚度为l的无限大平板，初温0度，初始瞬间将其放于温度为100度的流体中，流体与板面间的表面传热系数为一常数。

试确定在非稳态过程中板内的温度分布。

(1)显式差分法：图3 显式差分法(2)隐式差分法：图4 隐式差分法小结：显式格式仅当时格式是稳定的。

（其中称为网格比）隐式格式从k层求k+1层时，需要求解一个阶方程组。

而且隐式格式的稳定性对网格比没有要求，即为绝对稳定的。

三、Possion方程：取f=1,R=1图5差分法图6 误差小结：观察误差曲面，其绝对误差数量级为附Matlab程序：第1题：%===========================Boundary Value Problem 1clear;clc;A=[-2.01 1 0 0 0 0 0 0 0;1 -2.01 1 0 0 0 0 0 0;0 1 -2.01 1 0 0 0 0 0;0 0 1 -2.01 1 0 0 0 0;0 0 0 1 -2.01 1 0 0 0;0 0 0 0 1 -2.01 1 0 0;0 0 0 0 0 1 -2.01 1 0;0 0 0 0 0 0 1 -2.01 1;0 0 0 0 0 0 0 1 -2.01;];c1=[0.1;0.2;0.3;0.4;0.5;0.6;0.7;0.8;0.9];C=0.01*c1-1*[0;0;0;0;0;0;0;0;1];y=A\C;x=0:0.1:1;yn=[0;y;1];ye=2*(exp(x)-exp(-x))/(exp(1)-exp(-1))-x;figure(1);plot(x,yn,'*',x,ye);legend('numerical solution','exact solution')xlabel('x','fontsize',20);ylabel('y','fontsize',20);set(gca,'fontsize',18);figure(2);err=abs(ye'-yn);plot(x,err);legend('error')xlabel('x','fontsize',20);ylabel('y','fontsize',20);set(gca,'fontsize',18);第2题：%========================Boundary Value Problem 1_Explicit %显式clear;clcl=20;%板厚h=1;%步长J=l/h;T=50;%时间tao=2.5;%步长N=T/tao;%下面组合A矩阵a=0.2;lamda=tao/(h^2);zhu=1-2*a*lamda;ce=a*lamda;a00=ones(1,J-1);a0=diag(a00);A0=zhu*a0;a01=ones(1,J-2);a1=diag(a01,1);A1=ce*a1;a2=diag(a01,-1);A2=ce*a2;A=A0+A1+A2;u(:,1)=0; %板的初始温度for i=2:N+1u(1,i)=100-100*exp(-(i-1)*tao); %边界条件u(J+1,i)=100-100*exp(-(i-1)*tao); %边界条件end% g01=u(:,1);% g02=u(:,J);for k=1:N% g01=ce*g1(1,k);% g02=ce*g2(1,k);oo=zeros(J-3,1);g(:,k)=[ce*u(1,k);oo;ce*u(J+1,k)];u(2:end-1,k+1)=A*u(2:end-1,k)+g(:,k);endt=0:h:l;x=0:tao:T;mesh(x,t,u)xlabel('t','fontsize',20);ylabel('x','fontsize',20);zlabel('T','fontsize',20);set(gca,'fontsize',18);%========================Boundary Value Problem 1_2Implicit %隐式clear;clcl=20;%板厚h=1;%步长J=l/h;T=50;%时间tao=2.5;%步长N=T/tao;%下面组合A矩阵a=0.2;lamda=tao/(h^2);zhu=1+2*a*lamda;ce=-a*lamda;a00=ones(1,J-1);a0=diag(a00);A0=zhu*a0;a01=ones(1,J-2);a1=diag(a01,1);A1=ce*a1;a2=diag(a01,-1);A2=ce*a2;A=A0+A1+A2;u(:,1)=0; %板的初始温度for i=2:N+1u(1,i)=100-100*exp(-(i-1)*tao); %边界条件u(J+1,i)=100-100*exp(-(i-1)*tao); %边界条件end% g01=u(:,1);% g02=u(:,J);for k=1:N% g01=ce*g1(1,k);% g02=ce*g2(1,k);oo=zeros(J-3,1);g(:,k)=[ce*u(1,k);oo;ce*u(J+1,k)];u(2:end-1,k+1)=inv(A)*u(2:end-1,k)-inv(A)*g(:,k);endt=0:h:l;x=0:tao:T;mesh(x,t,u)xlabel('t','fontsize',20);ylabel('x','fontsize',20);zlabel('T','fontsize',20);set(gca,'fontsize',18);第3题：%=============================used by centered difference clear;function pdemodel[pde_fig,ax]=pdeinit;pdetool('appl_cb',1);set(ax,'DataAspectRatio',[1 1 1]);set(ax,'PlotBoxAspectRatio',[1.5 1 1]);set(ax,'XLim',[-1.5 1.5]);set(ax,'YLim',[-1 1]);set(ax,'XTickMode','auto');set(ax,'YTickMode','auto');% Geometry description:pdecirc(0,0,1,'C1');set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','C1')% Boundary conditions:pdetool('changemode',0)pdesetbd(4,...'dir',...1,...'1',...'0')pdesetbd(3,...'dir',...1,...'1',...'0')pdesetbd(2,...'dir',...1,...'1',...'0')pdesetbd(1,...'dir',...1,...'1',...'0')% Mesh generation:setappdata(pde_fig,'Hgrad',1.3);setappdata(pde_fig,'refinemethod','regular');pdetool('initmesh')pdetool('refine')% PDE coefficients:pdeseteq(1,...'1.0',...'0.0',...'1',...'1.0',...'0:10',...'0.0',...'0.0',...'[0 100]')setappdata(pde_fig,'currparam',...['1.0';...'0.0';...'1 ';...'1.0'])% Solve parameters:setappdata(pde_fig,'solveparam',...str2mat('0','1524','10','pdeadworst',...'0.5','longest','0','1E-4','','fixed','Inf'))% Plotflags and user data strings:setappdata(pde_fig,'plotflags',[1 1 1 1 1 1 1 1 0 0 0 1 0 0 1 0 0 1]); setappdata(pde_fig,'colstring','');setappdata(pde_fig,'arrowstring','');setappdata(pde_fig,'deformstring','');setappdata(pde_fig,'heightstring','');% Solve PDE:pdetool('solve')。

高维波动方程数值模拟的隐式分裂有限差分格式

（ｚ）一ｐ（ｘ＋Ａｘ）一２ｐ（ｘ）＋ｐ（ｘ — Ａｘ）
辅助变量Ｕ，Ｖ和Ｗ定义为
ｆＵ一嵋一２ＰＴｉｋ＋ｕ％ｏ
这就要求在实际计算中选取较小的模拟时间步长
以确保数值稳定性，从而降低了计算效率。我们提出一种新的、与文献［１］类似的隐式分裂有限差分
高维波动方程数值模拟的隐式分裂有限差分格式597图5复杂的盐丘模型图6八阶显式差分格式的瞬时波场图7隐式分裂有限差分格式的瞬时波场结束语本文隐式分裂有限差分格式求解多维波动方程的基本原理是将多维波动方程分解为一维方程用追赶法快速地求解
维普资讯
第４６卷第６期２００７年１１月
格式，对于给定的高维波动方程，首先按照不同的
传播方向进行分解，然后对分解后的一维方程分别
｛一硪Ｌ２ｐＴｊ￣＋ｌ一磁一２ｐＴｊｋ＋
更新后的波场户可由
（３）
硅一盟
的一维波动问题，然后分别沿各方向隐式求解。该格式包含了，Ｙ，ｚ三个方向相互独立的一维隐式差分格式，每个方向的一维格式在数值离散后归结为一个三对角矩阵问题，可以用追赶法快速地求解。将该格式从时间一
空间域变换至时间一波数域，证明此格式可以通过适当地选取参数来提高计算精度，保证计算过程的稳定性和与
式中：ｐ（ｘ，Ｙ，ｚ）是压力波场；ｆ（ｚ，Ｙ，ｚ）是已知的速度场。为了对（１）式进行差分求解，将ｐ（ｘ，Ｙ，ｚ）离

基于交变隐式差分方向方法的时域有限差分法

基于交变隐式差分方向方法的时域有限差分法摘要本文主要针对基于交变隐式差分方向方法的时域有限差分法（Alternating Direction Implicit Finite Difference Time Domain method，简称ADI-FDTD方法）做了一定研究。

论文首先介绍了二维ADI-FDTD方法，就其数值稳定性和数值色散特性进行了研究，验算了ADI-FDTD的Mur吸收边界条件对其稳定性的影响。

关键词：基于交变隐式差分方向方法的时域有限差分法有限时域差分法数值稳定数值色散吸收边界条件一、二维ADI-FDTD 方法基本原理基于交变隐式差分方向方法的时域有限差分法也就是ADI-FDTD 方法。

传统的FDTD 方法,属于显示差分方法,因此具有显示差分方法的共同特性。

其离散格式有以下两个方面的特点：一个就是数值色散对空间离散网格的要求,空间离散网格尺寸必须为所要模拟电磁波最短波长的1/12，通常在程序中取1/20以减少数值计算带来的误差。

第二就是时间和空间离散间隔之间应当满足Courant-Friedrich-Levy(CFL)稳定性条件，或者简称Courant 稳定条件，即：222max )(1)(1)(11z y x t v ∆+∆+∆≤∆其中Vmax 为电磁波在媒质中传播的最大相速。

如果时间步长不满足上述条件，将导致FDTD 的计算发散，特别是目标较之入射波的波长有细微结构时，随着空间网格尺寸变小，为了满足稳定性条件，时间步长也相应的取小，致使总的CPU 计算时间有可能达到无法实现的地步。

为了提高FDTD 方法的计算效率和应用范围，90年代后期提出了多种与其他技术相结合的混合方法。

第六章就介绍最后一种基于交变隐式差分方向方法的时域有限差分法，也就是ADI-FDTD 方法。

与显式差分方法相反，隐式差分格式总是稳定的，其时间步长仅受数值误差的限制。

但是，隐式差分格式的缺点是需要通过矩阵求逆，或者是迭代求解大型线性方程组，计算复杂且量大。

水环境保护作业代码

水环境保护课程作业——隐式差分法计算河段BOD浓度【例题】某均匀河段长8 km，流速u =5km/h，纵向离散系数E =2km2/h，BOD的降解系数K1=0.0151h-1，上游断面有一断面混合均匀的污染源，0~1h间稳定排污，使该断面的污染浓度保持在20mg/L，以后排污停止，试用隐式差分法计算本河段各断面的BOD浓度变化过程。

【解题思路】1、建立所需变量和数组。

经分析变量分别有纵向离散系数E、流速u、BOD的降解系数K、河段长L、时间步长△t和距离步长△d，用于进行循环的i，j。

取时间步长为0.1h，距离步长为0.5km，可得需要一11行17列的数组来存放所得数据结果，即不同时间点不同位置的浓度值。

2、进行边界的初值的赋值：第一行和第一列的数值为已知条件3、利用追赶法顺序求解不同时刻的g和w的值，根据公式的不同，分两步求解，第一步是g[1]和w[1]的值，第二步是循环求解第二个到倒数第二个的g和w的值，每个时刻的最后的一个g值等于最后的浓度值。

4、利用追赶法倒序求出每个时刻，各个位置点的浓度值。

思路框图如下：【编程代码】#include <stdio.h>void main(){double k=0.0151,d=0.5;int u=5,e=2;int i,j,m=10,n=16;double a=0,b=0,r=0,t=0.1;doublex[17]={0},g[17]={0},w[17]={0},aa[11][ 17]={0};//基础变量和数组的定义FILE *fp;if((fp=fopen("D:\\水环境保护.txt","w"))==NULL){printf("can not open file\n");return;}//结果的输出a=-e/(d*d);b=1/t+2*e/(d*d)+k/2;r=a;for(i=1;i<=n;i++){aa[0][i]=0;}for(i=0;i<m;i++){aa[i][0]=20;}aa[m][0]=0;//对边界初值的赋值for(j=0;j<=m;j++){w[1]=r/b;for(i=1;i<=n;i++){x[i]=aa[j][i]*(1/t-u/d)+aa[j][i-1]*(u/ d-k/2);}x[1]-=a*aa[j+1][0];g[1]=x[1]/b;for(i=2;i<n;i++){g[i]=(x[i]-a*g[i-1])/(b-a*w[i-1]);w[i]=r/(b-a*w[i-1]);}//求w[i]g[i]的数值g[16]=x[16]/(b+2*r);aa[j+1][n]=g[n];for(i=15;i>0;i--){aa[j+1][i]=g[i]-w[i]*zl[j+1][i+1];}//求各点的浓度值for(i=0;i<=n;i++){printf("%7.3f",aa[j][i]);fprintf(fp,"%7.3f",aa[j][i]);}printf("\n");}}【数据结果】【图表】。

【毕业设计(论文)】二维热传导方程有限差分法的MATLAB实现

第1章前言1.1问题背景在史策教授的《一维热传导方程有限差分法的MATLAB实现》和曹刚教授的《一维偏微分方程的基本解》中，对偏微分方程的解得MATLAB实现问题进行过研究，但只停留在一维中，而实际中二维和三维的应用更加广泛。

诸如粒子扩散或神经细胞的动作电位。

也可以作为某些金融现象的模型，诸如布莱克-斯科尔斯模型与Ornstein-uhlenbeck过程。

热方程及其非线性的推广形式也被应用与影响分析。

在科学和技术发展过程中，科学的理论和科学的实验一直是两种重要的科学方法和手段。

虽然这两种科学方法都有十分重要的作用，但是一些研究对象往往由于他们的特性（例如太大或太小，太快或太慢）不能精确的用理论描述或用实验手段来实现。

自从计算机出现和发展以来，模拟那些不容易观察到的现象，得到实际应用所需要的数值结果，解释各种现象的规律和基本性质。

科学计算在各门自然科学和技术科学与工程科学中其越来越大的作用，在很多重要领域中成为不可缺少的重要工具。

而科学与工程计算中最重要的内容就是求解科学研究和工程技术中出现的各种各样的偏微分方程或方程组。

解偏微分方程已经成为科学与工程计算的核心内容，包括一些大型的计算和很多已经成为常规的计算。

为什么它在当代能发挥这样大的作用呢？第一是计算机本身有了很大的发展；第二是数值求解方程的计算法有了很大的发展，这两者对人们计算能力的发展都是十分重要的。

1.2问题现状近三十年来，解偏微分方程的理论和方法有了很大的发展，而且在各个学科技术的领域中应用也愈来愈广泛，在我国，偏微分方程数值解法作为一门课程，不但在计算数学专业，而且也在其他理工科专业的研究生的大学生中开设。

同时，求解热传导方程的数值算法也取得巨大进展，特别是有限差分法方面，此算法的特点是在内边界处设计不同于整体的格式，将全局的隐式计算化为局部的分段隐式计算。

而且精度上更好。

目前，在欧美各国MATLAB的使用十分普及。

在大学的数学、工程和科学系科，MATLAB苏佳园：二维热传导方程有限差分法的MATLAB实现被用作许多课程的辅助教学手段，MATLAB也成为大学生们必不可少的计算工具，甚至是一项必须掌握的基本技能。

有限差分编程书籍

有限差分编程书籍
有限差分编程是一种常用的数值计算方法，以下是一些关于有限差分编程的书籍推荐：- 《Computational fluid dynamics, principles and applications》：由北大工学院的几位老教授编写，该书对经典的格式讲解得非常清楚，具有很强的实用性。

- 《Computational Fluid Dynamics, the basics with applications》：此书被很多人推荐，虽然其中的算法基本已经过时，但作为基础教材学习还是不错的。

- 《The finite-volume method method in computational fluid dynamics, an advanced introduction to Openfoam and Matlab》：这本书是目前有限体积编程方面的优质教材，其优势包括：①所讲算法具有实用性，且全部为非结构网格算法；②有大量C语言的代码；
③图片示意清晰，方便初学者理解。

有限差分编程的书籍还有很多，你可以根据自己的需求和水平选择适合的书籍进行学习。

隐式差分方程课件

1 2
k 2 Dx4

)u
n 1 m

u
n m
式中左边如果仅保留二阶导数项，且以
格式
(1
k h2

2 x
)U
n1 m

U
n m
1 h2

2 x
替代 Dx2
，则得差分
或者
rU
n1 m1

(1
2r
)U
n1 m

rU
n1 m1

U
n m
（2.41）
格式用图2.5表示，其截断误差阶为 (k2 h2) ,与古典差分格式相同。
exp(
1 2
k L)umn1

exp(1 2
k L)umn
由 L Dx2
得 [1
1 2
k Dx2

1 2
(1 2
kDx2 )2
]umn1

[1
1 2
k Dx2

1 2
(1 2
kDx2 )2
]umn
（2.42）
两边仅保留前二项，用
1 h2

代替 2
x
Dx2
，则得差分格式
(1
232cranknicolson隐式格式cranknicolson隐式差分格式是解热传导方程226的常用的差分格式为了推导它由式224有242两边仅保留前二项用代替则得差分格式243这是一个隐式差分格式称为cranknicolson差分格式截断误差阶为也可写为kdkdkdkd244由于格式244中包括六个结点故也可称为六点格式如图26所示
0.956 821 703 419 0.000 079

有限差分法

有限差分法有限差分法finite difference method微分方程和积分微分方程数值解的方法。

基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。

然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。

有限差分法的主要内容包括：如何根据问题的特点将定解区域作网格剖分；如何把原微分方程离散化为差分方程组以及如何解此代数方程组。

此外为了保证计算过程的可行和计算结果的正确，还需从理论上分析差分方程组的性态，包括解的唯一性、存在性和差分格式的相容性、收敛性和稳定性。

对于一个微分方程建立的各种差分格式，为了有实用意义，一个基本要求是它们能够任意逼近微分方程，这就是相容性要求。

另外，一个差分格式是否有用，最终要看差分方程的精确解能否任意逼近微分方程的解，这就是收敛性的概念。

此外，还有一个重要的概念必须考虑，即差分格式的稳定性。

因为差分格式的计算过程是逐层推进的，在计算第n＋1层的近似值时要用到第n层的近似值，直到与初始值有关。

前面各层若有舍入误差，必然影响到后面各层的值，如果误差的影响越来越大，以致差分格式的精确解的面貌完全被掩盖，这种格式是不稳定的，相反如果误差的传播是可以控制的，就认为格式是稳定的。

只有在这种情形，差分格式在实际计算中的近似解才可能任意逼近差分方程的精确解。

关于差分格式的构造一般有以下3种方法。

最常用的方法是数值微分法，比如用差商代替微商等。

另一方法叫积分插值法，因为在实际问题中得出的微分方程常常反映物理上的某种守恒原理，一般可以通过积分形式来表示。

此外还可以用待定系数法构造一些精度较高的差分格式。