扩散方程的差分解法

扩散方程的差分解法

在研究热传导过程、扩散过程、边界层现象时，我们常常遇到抛物型方程，这类方程中最典型、最简单的就是热传导方程。热传导方程中的自变量中包括时间t ，它是描述一种随时间变化的物理过程，即所谓不定常现象。这类问题的基本定解问题应是初值问题，即在初始时刻（t=0）时给定定解条件，求解t>0时的解。

本文主要运用有限差分法对一维扩散方程进行求解，并对差分解的适定性、相容性、收敛性及稳定性进行分析，同时与解析解进行对比。

1.扩散方程

一维扩散方程为：

22u u

t x

α??=??

（1）

式中，u 为因知量，α为扩散系数，x 为坐标，t 为时间。

其定解条件如下：

初始条件：

(,0)() 0x u x f x L

=≤≤

（2）

边界条件：

12(0,)() , (,)()u t f t u L t f t ==

（3） 一般假定函数()f x ，1()f t ，2()f t 满足连接条件，即1(0)(0) f f =，2()(0) f L f =。

2.有限差分法

有限差分法是数值计算解微分方程古老的方法之一，也是系统化地、数值地求解数学物理方法的方程。其控制方程中的导数用离散点上函数值的差商代替。

差分格式可以分为显格式和隐格式。所谓显格式是指在任一结点上因变量在新是时间层上的值可以通过之前的时间层上相邻结点变量的值显式解出来。由于这些层的变量值是已知的，当时间向前推进时，空间点上的新的变量值就只需逐点计算就行了，因此显格式计算起来比较省事。隐格式则是指任一结点上变量在新的时间层的值，不能通过之前的时间层上相邻结点的值显式解出来，它不仅与之前的时间层上的已知值有关，而且也与新时间层的相邻结点的变量值有关。因而一个差分方程常常包括几个相邻结点上的未知数，未知数的个数取决于格式的构成形式。为了解出这些未知数需要联立新的方程，而每引进一个新的方程往往又同时引进了新的未知数。因此，隐格式总是伴随着求解巨大的代数方程组。隐格式的主要缺点是计算工作量大，因而不如显格式计算得快，但这只是就时间步长一样的情况而言的。隐格式的主要优点是时间步长可以比显格式能够采用的最大步长大很多。显格式的时间步长受到稳定性条件的限制，而隐格式则几乎不受限制。

3.方程的离散

3.1 显格式

采用时间前差及第n 时间层的空间中心差，得一维扩散方程的显格式解：

111

2

2()n n n n n

j j

j j j u u u u u t

x α

++---+=??

（4）

即

111(2)

n n n n n

j j j j j u u r u u u ++-=+-+

（5）

式中，2()t r x α

?=?

3.2 全隐格式

采用时间前差及第（n+1）时间层的空间中心差，得一维扩散方程的全隐格式解：

1111

11

2

2()n n n n n j j

j j j u u u u u t

x α

+++++---+=??

（6）

4.差分解的基本问题

差分解的基本问题包括：适定性、相容性、收敛性和稳定性四个方面。

4.1 适定性

在用差分方程作微分方程数值解时，首先，要求微分方程的问题是适定的。所谓适定性问题是指这一微分方程在一定的初始和边界条件下要有唯一解，并且在初始条件和边界条件稍有改变时，微分方程的解也只是稍有偏离。从数学的角度讲，若微分方程的解存在并且是唯一的，同时连续依赖于数据（初始条件、边界条件），则问题是适定的。

对于抛物线方程，这是初值问题，要求给出初始条件，并且在区域边界上给出边界条件。在本文的一维扩散问题中，除了给出t=0时的初值(,0)u x 外，应在左、右两端边界，即0x =及x L =，各给定一个边界条件，第一类边界条件u 的值，或给定第二类边界条件u n ?的值或给定第三类边界条件u 与u n ??的组合。

由于抛物型方程具有扩散性质，它往往是随时间扩散的，u 值将不断因扩散而衰减。 在本文的一维扩散方程中，给定了初始条件，同时在区域的左、右端边界给定了边界条件，满足适定性要求。

4.2 相容性

将一个偏微分方程用差分格式化为相应差分方程，当步长t ?和x ?趋近于零时，这个差分方程应当收敛于原微分方程，也就是说，相应的差分方程和微分方程之间的截断误差在任一时刻任一网格点上均应趋近于零，这样的差分方程和微分方程才是相容的。 对一维扩散方程：

220

u u Lu t x α??=-=??

（7）

采用显格式差分格式，令t τ?=，x h ?=，则

1112

20

n n n n n j j j j j n

h j u u u u u L u h ατ++---+=-= （8）

用Taylor 级数展开代入上述差分方程中，则有

()222424

224(),0

212n

n

n n h j j j j u u u h u L u O h t x t x ταατ????????=-+-+= ? ?????????

（9）

当0τ→，0h →时，上式化为

22()0

n

n

h j

j u u L u t x α??=-=??

（10） 可见，此差分格式所构成的差分方程与原来的微分方程是相容的，故该显格式为相容格式。

采用全隐格式差分格式，令t τ?=，x h ?=，则

111111

2

20n n n n n j j

j j j n h j

u u u u u L u h α

τ

+++++---+=

-=

（11）

用Taylor 级数展开代入上述差分方程中，则有

()222(),0

n

n

h j

j u u L u O h t x ατ??=-+=??

（12）

当0τ→，0h →时，上式化为

22()0

n n

h j

j u u L u t x α??=-=??

（13） 可见，此差分格式所构成的差分方程与原来的微分方程是相容的，故该全隐格式为相容格式。

4.3 稳定性

差分法计算中所产生的误差（舍入误差，参数误差等）随时间衰减或不增大，则称离散格式是稳定的，反之，则是不稳定的。分析差分方程稳定性有不同的方法，如矩阵方法，谐波分析法等。下面用谐波分析法对以上两种离散格式的稳定性进行分析： （1）显格式 显格式的差分方程为

111

2

20

()n n n n n

j j

j j j n h j

u u u u u L u t

x α

++---+=

-=??

（14）

即

111(12)()

n n n n j j j j u r u r u u ++-=-++

（15）

其误差方程为

111(12)()

n n n n j j j j r r εεεε++-=-++

（16）

任取一k 次谐波分量 e

e j

ikx n ikjh

j k k A A ε==

（17）

则

()1e (12)e e (122cos )

n ikjh ikh ikh n

j k j A r r r r kh εε+-??==-++=-+??

（18）

误差放大因子为

122cos 12(1cos )

r r kh r kh λ=-+=--

（19） 要满足稳定性条件，则要求对所有的k 值均有||1λ≤，须|14|1r -≤，即1

02r ≤≤。

因此，一维扩散方程显格式的稳定性条件为

2

(),0

2x t αα

??≤≥且 （20）

（2）全隐格式

全隐格式的差分方程为

1111

11

2

20()n n n n n j j

j j j n h j

u u u u u L u t

x α

+++++---+=

-=??

（21）

即

111

11(12)()

n n n n j j j j r u u r u u ++++-+=++

（22）

其误差方程为

111

11(12)()

n n n n j j j j r r εεεε++++-+=++

（23）

任取一k 次谐波分量 1

e e j

ikx n ikjh

j

k k A A ε+==

（24）

则

111e n n ikh

j j

εε+++=，

111e n n ikh

j j

εε++--=

（25）

则误差方程为

11

[(12)r(e e )][12(1cos )]ikh ikh n n n j j j

r r kh εεε-+++-+=+-=

（26）

误差放大因子为

1

12(1cos )]

r kh λ=

+-

（27）

要满足稳定性条件，则要求对所有的k 值均有||1λ≤。从（28）式中可以看出，当0r ≥（即0α≥）时，||1λ≤恒成立。因此，全隐格式是无条件稳定的。

4.4 收敛性

如果差分方程的解为j u ，微分方程的解为u ，若当0t ?→，0x ?→时，差分方程的解与微分方程的解之差

()0

j j u u x ε=-→ （28）

则称差分格式是收敛的。

拉克斯（Lax ）等价定理指出，如果问题是适定的，并且差分格式满足相容性条件，那么差分格式的稳定性就是该格式收敛性的充分而必要的条件。

由4.1、4.2和4.3的分析可知，显格式和全隐格式都满足适定性、相容性及稳定性的条件，因而这两种格式满足收敛性要求。

5.差分方程的求解

5.1 初始条件及边界条件

由以上对一维扩散问题的分析，可知，求解一维扩散方程需给定初始条件及边界条件。 在本文计算中，取10L m =，1α=。 初始条件（0t =时）

2/ 0/2

(,0)()

2(1/) /2x L x L u x f x x L L x L ≤≤?==?-≤≤?

（29）

边界条件为

12(0,)()0

0(,)()0

u t f t t u L t f t ==?>?==? （30）

其初始时刻（0t =）时的u 分布如图1所示，x=0m 处u 随时间变化情况如图2所示，x=10m

处u 随时间变化情况如图3所示。

图1 初始时刻u分布图

图2 x=0处u随时间变化图

图3 x=10m处u随时间变化图

5.2 空间步长及时间步长

通过对差分格式的稳定性分析知，显格式空间步长与时间步长间应满足一定的关系，即式（21）。本文选用的步长如表1所示，分别用显格式和全隐格式进行数值计算，差分网格如图4所示。

图4 差分网格示意图

5.2 求解方法

对于显格式，有111(2)n n n n n

j j j j j u u r u u u ++-=+-+，若n 时间层上的u 已知，则可直接求解出

（n+1）时间层上的u 值。因此，由于给出了初始条件，显格式无需迭代，直接从t=0时刻

开始，逐一计算下一个时间层上的u 值，便可求解出各时间层上各空间点的u 值。

对于全隐格式，由于1n j u +与11n j u ++，1

1n j u +-有联系，不能直接求解，必须联解代数方程组。

此时方程组为三对角方程组，可采用追赶法进行求解。 5.3 计算程序 5.3.1显格式

!----------------------------------------显格式求解扩散方程---------------------------------------------------! !----------------------参数设置-----------------------------------------！ parameter (nt=20001,nx=101) !nt:时间节点数;nx:空间节点数 dimension u(nx,nt)

Sx=10.0!计算长度（m ） dt=0.003!时间步长（s ） dx=0.1!空间步长（m ） alfa=1.0!扩散系数

!-----------!

open(1,file='显格式沿程变化.txt')

open(2,file='显格式随时间变化.txt')

r=alfa*dt/dx**2

!----------------------边界条件----------------------------------------!

do n=1,nt

u(1,n)=0

u(nx,n)=0

enddo

!----------!

!---------------------初始条件-----------------------------------------!

do j=1,nx

x=(j-1)*dx

if (x.le.sx/2) then

u(j,1)=2*x /sx

else

u(j,1)=2*(1-x/sx)

endif

enddo

!-----------!

!--------------------求解差分方程------------------------------------!

do n=2,nt

do j=2,nx-1

u(j,n)=u(j,n-1)+r*(u(j-1,n-1)-2*u(j,n-1)+u(j+1,n-1)) enddo

enddo

!-----------!

!-------------------每隔15s输出沿程变化------------------------!

do n=1,nt

t=(n-1)*dt

if (mod(int(t*1000),15000).eq.0) then

write(1,*) 't=',t,'s'

do j=1,nx

write(1,*) (j-1)*dx,u(j,n)

enddo

endif

enddo

!-----------!

!-------------------每隔2.5m输出随时间变化------------------------! do j=1,nx

x=dx*(j-1)

if (mod(int(x*10),25).eq.0) then

write(2,*) 'x=',x,'m'

do n=1,nt,200

write(2,*) (n-1)*dt,u(j,n)

enddo

endif

enddo

!-----------!

end

5.3.2全隐格式

!----------------------------------------全隐格式求解扩散方程-----------------------------------------------! !----------------------参数设置-----------------------------------------！

!implicit double precision (a-h,o-z)

parameter (nt=20001,nx=101) !nt:时间节点数;nx:空间节点数

dimension u(nx,nt),a(nx),b(nx),c(nx-1),d(nx-1),xx(nx)

Sx=10!计算长度（m）

dt=0.003!时间步长（s）

dx=0.1!空间步长（m）

alfa=1.0!扩散系数

!-----------!

open(1,file='全隐格式沿程变化.txt')

open(2,file='全隐格式随时间变化.txt')

r=alfa*dt/dx**2

!----------------------边界条件----------------------------------------!

do n=1,nt

u(1,n)=0

u(nx,n)=0

enddo

!----------!

!---------------------初始条件-----------------------------------------!

do j=1,nx

x=(j-1)*dx

if (x.le.sx/2) then

u(j,1)=2*x/sx

else

u(j,1)=2*(1-x/sx)

endif

enddo

!-----------!

!-----------------------------三对角方程组系数矩阵-------------------------! a(1)=1

do i=1,nx-1

d(i)=-r

c(i)=-r

a(i+1)=1+2*r

enddo

c(1)=0

!-----------!

!-----------------------------三对角方程组常数项-------------------------! do i=1,nx

b(i)=u(i,1)

enddo

!-----------!

!----------------------------------求解差分方程-----------------------------! do n=2,nt

call systri(nx,d,a,c,b,xx)

do j=1,nx-1

b(j)=xx(j)

u(j,n)=xx(j)

enddo

enddo

!-----------!

!-------------------每隔15s输出沿程变化------------------------!

do n=1,nt

t=(n-1)*dt

if (mod(int(t*1000),15000).eq.0) then

write(1,*) 't=',t,'s'

do j=1,nx

write(1,*) (j-1)*dx,u(j,n)

enddo

endif

enddo

!-----------!

!-------------------每隔2.5m输出随时间变化------------------------!

do j=1,nx

x=dx*(j-1)

if (mod(int(x*10),25).eq.0) then

write(2,*) 'x=',x,'m'

do n=1,nt,200

write(2,*) (n-1)*dt,u(j,n)

enddo

endif

enddo

!-----------!

end

subroutine systri(n,d,a,c,b,x)

! d：下对角线

! c：上对角线

! a：主对角线

! b：方程右边常数项

! x：方程的解

! implicit double precision (a-h,o-z)

dimension a(n),d(n-1),c(n-1)

dimension p(n),q(n-1),y(n),x(n),b(n)

p(1)=a(1)

do i=1,n-1

q(i)=c(i)/p(i)

p(i+1)=a(i+1)-d(i)*q(i)

enddo

y(1)=b(1)/p(1)

do i=2,n

y(i)=(b(i)-d(i-1)*y(i-1))/p(i)

enddo

x(n)=y(n)

do i=n-1,1,-1

x(i)=y(i)-q(i)*x(i+1)

enddo

return

end

6 计算结果与分析

6.1显格式计算结果与分析

6.1.1随时间变化

经计算，沿程断面u随时间变化如下图5-图7所示，从图中可以看出，各断面u值随着时间逐渐递减，在t=60s时，u值衰减至接近于0，曲线的斜率随着时间逐渐趋缓，可见u 的衰减速率逐渐减小。由初始条件可知，初始时刻，u分布关于x=5m对称，对比x=2.5m及x=7.5m的u值随时间变化可知，这两个断面的u随时间变化关系一致，说明对称分布的初始条件随着时间的变化，其分布仍具有对称性。

图5 x=2.5m处u随时间变化

图6 x=5m处u随时间变化

图7 x=7.5m处u随时间变化

6.1.2沿程分布

图8-图11给出了t=15s、t=30s、t=45s及t=60s时u值的分布情况，从图中可以看出，u 值的分布呈抛物线分布，在x=5m处达到最大值，同时u分布关于x=5m对称。可见，u分布仍保持初始时刻的对称关系，只是分布曲线趋于平缓。对比各时刻u的峰值，可知，u的峰值随着时间逐渐减小，至t=60s时，其峰值为0.002左右，此时最大值与最小值之差在0.002左右，可见，在t=60s时，扩散趋于稳定，即整个区域u值基本相同。

图8 t=15s时u沿程分布

图9 t=30s时u沿程分布

图10 t=45s时u沿程分布

图11 t=60s时u沿程分布

6.2全隐格式计算结果与分析

6.2.1随时间变化

经计算，沿程断面u随时间变化如下图12-图14所示，从图中可以看出，各断面u值随着时间逐渐递减，在t=60s时，u值衰减至接近于0，曲线的斜率随着时间逐渐趋缓，可见u的衰减速率逐渐减小。由初始条件可知，初始时刻，u分布关于x=5m对称，对比x=2.5m 及x=7.5m的u值随时间变化可知，这两个断面的u随时间变化关系一致，说明对称分布的初始条件随着时间的变化，其分布仍具有对称性。

图12 x=2.5m处u随时间变化

图13 x=5.0m处u随时间变化

图14 x=7.5m处u随时间变化

6.2.2沿程分布

图15-图18给出了t=15s、t=30s、t=45s及t=60s时u值的分布情况，从图中可以看出，u值的分布近似呈抛物线分布，在x=5m处达到最大值，同时u分布关于x=5m对称。可见，u分布仍保持初始时刻的对称关系，只是分布曲线趋于平缓。对比各时刻u的峰值，可知，u的峰值随着时间逐渐减小，至t=60s时，其峰值为0.002左右，此时最大值与最小值之差在0.002左右，可见，在t=60s时，扩散趋于稳定，即整个区域u值基本相同。

图15 t=15s时u沿程分布

图16 t=30s时u沿程分布

图17 t=45s时u沿程分布

图18 t=60s时u沿程分布

6.3计算结果对比

6.3.1随时间变化对比

图19-图21给出了x=2.5m，x=5m，x=7.5m处显格式与全隐格式计算所得u值随时间变化的对比图，从图中可以看出，显格式与全隐格式计算得到的曲线几乎重合基本一致，说明两者计算结果基本一致。两种格式计算所得的各断面u值均随着时间逐渐递减，在t=60s 时，u值均衰减至接近于0，曲线的斜率均随着时间逐渐趋缓，可见两种格式计算中u的衰减速率均逐渐减小。由初始条件可知，初始时刻，u值分布关于x=5m对称，对比x=2.5m及x=7.5m的u值随时间变化可知，在两种计算格式下，这两个断面的u值随时间变化关系一致，说明对称分布的初始条件在两种计算格式下，均能很好地保持对称性。

图19 x=2.5m处随时间变化对比

图20 x=5m处随时间变化对比

图21 x=7.5m处随时间变化对比

6.3.2沿程分布对比

图22-图25给出了t=15s、t=30s、t=45s及t=60s时u值的沿程分布对比图，从图中可以看出，两种计算格式下u值的分布均近似呈抛物线分布，在x=5m处达到最大值，同时u值分布均关于x=5m对称。可见，两种计算格式得到的u值分布仍保持初始时刻的对称关系，只是分布曲线趋于平缓。对比各时刻u的峰值，可知，u的峰值随着时间逐渐减小，至t=60s 时，其峰值为0.002左右，此时最大值与最小值之差在0.002左右，可见，在t=60s时，扩散趋于稳定，即整个区域u值基本相同。另一方面，从图中可以看出全隐格式计算得到的u 值较显格式的略大，说明，全隐格式计算的扩散较显格式计算的为慢。

图22 t=15s时沿程分布对比

图23 t=30s时沿程分布对比

图24 t=45s 时沿程分布对比

图25 t=60s 时沿程分布对比

6.1.3对比解析解

对于本文研究问题，其有分析解为

()22221811sin sin exp()2n u n n n t n ππππ∞=??

=-∑ ???

（31）

对比x=5m 处的分析解与显格式差分解如表2所示。

表2 分析解与差分解对比

(完整版)大连理工大学高等数值分析抛物型方程有限差分法

抛物型方程有限差分法 1. 简单差分法 考虑一维模型热传导方程 (1.1) )(22x f x u a t u +??=??，T t ≤<0 其中a 为常数。)(x f 是给定的连续函数。（1.1）的定解问题分两类： 第一，初值问题（Cauchy 问题）：求足够光滑的函数()t x u ,，满足方程（1.1）和初始条件： （1.2） ()()x x u ?=0,, ∞<<∞-x 第二，初边值问题（也称混合问题）：求足够光滑的函数()t x u ,，满足方程（1.1）和初始条件： ()13.1 ()()x x u ?=0,, l x l <<- 及边值条件 ()23.1 ()()0,,0==t l u t u ， T t ≤≤0 假定()x f 和()x ?在相应的区域光滑，并且于()0,0，()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。

现在考虑边值问题（1.1），（1.3）的差分逼近 取 N l h = 为空间步长，M T = τ为时间步长，其中N ，M 是 自然数， jh x x j ==, ()N j ,,1,0Λ=； τ k y y k ==, ()M k ,,1,0Λ= 将矩形域G {}T t l x ≤≤≤≤=0;0分割成矩形网格。其中 ()j i y x ,表 示网格节点； h G 表示网格内点（位于开矩形G 中的网格节点）的集合； h G 表示位于闭矩形G 中的网格节点的集合； h Γ表示h G -h G 网格边界点的集合。 k j u 表示定义在网点()k i t x ,处的待求近似解，N j ≤≤0，M k ≤≤0。 注意到在节点()k i t x ,处的微商和差商之间的下列关系 （(,)k j k j u u x t t t ????≡ ? ????）： ()() ()ττ O t u t x u t x u k j k j k j +??? ????=-+,,1 ()() ()2112,,ττ O t u t x u t x u k j k j k j +??? ????=--+ ()()()h O x u h t x u t x u k j k j k j +??? ????=-+,,1 ()() ()h O x u h t x u t x u k j k j k j +??? ????=--,,1 ()() ()2112,,h O x u h t x u t x u k j k j k j +??? ????=--+ ()()() ()2 222 11,,2,h O x u h t x u t x u t x u k j k j k j k j +???? ????=+--+ 可得到以下几种最简差分格式

扩散方程的差分解法

扩散方程的差分解法 在研究热传导过程、扩散过程、边界层现象时，我们常常遇到抛物型方程，这类方程中最典型、最简单的就是热传导方程。热传导方程中的自变量中包括时间t ，它是描述一种随时间变化的物理过程，即所谓不定常现象。这类问题的基本定解问题应是初值问题，即在初始时刻（t=0）时给定定解条件，求解t>0时的解。 本文主要运用有限差分法对一维扩散方程进行求解，并对差分解的适定性、相容性、收敛性及稳定性进行分析，同时与解析解进行对比。 1.扩散方程 一维扩散方程为： 22u u t x α??=?? （1） 式中，u 为因知量，α为扩散系数，x 为坐标，t 为时间。 其定解条件如下： 初始条件： (,0)() 0x u x f x L =≤≤ （2） 边界条件： 12(0,)() , (,)()u t f t u L t f t == （3） 一般假定函数()f x ，1()f t ，2()f t 满足连接条件，即1(0)(0) f f =，2()(0) f L f =。 2.有限差分法 有限差分法是数值计算解微分方程古老的方法之一，也是系统化地、数值地求解数学物理方法的方程。其控制方程中的导数用离散点上函数值的差商代替。 差分格式可以分为显格式和隐格式。所谓显格式是指在任一结点上因变量在新是时间层上的值可以通过之前的时间层上相邻结点变量的值显式解出来。由于这些层的变量值是已知的，当时间向前推进时，空间点上的新的变量值就只需逐点计算就行了，因此显格式计算起来比较省事。隐格式则是指任一结点上变量在新的时间层的值，不能通过之前的时间层上相邻结点的值显式解出来，它不仅与之前的时间层上的已知值有关，而且也与新时间层的相邻结点的变量值有关。因而一个差分方程常常包括几个相邻结点上的未知数，未知数的个数取决于格式的构成形式。为了解出这些未知数需要联立新的方程，而每引进一个新的方程往往又同时引进了新的未知数。因此，隐格式总是伴随着求解巨大的代数方程组。隐格式的主要缺点是计算工作量大，因而不如显格式计算得快，但这只是就时间步长一样的情况而言的。隐格式的主要优点是时间步长可以比显格式能够采用的最大步长大很多。显格式的时间步长受到稳定性条件的限制，而隐格式则几乎不受限制。 3.方程的离散 3.1 显格式 采用时间前差及第n 时间层的空间中心差，得一维扩散方程的显格式解： 111 2 2()n n n n n j j j j j u u u u u t x α ++---+=?? （4） 即 111(2) n n n n n j j j j j u u r u u u ++-=+-+ （5）

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

第1章前言 1.1问题背景 在史策教授的《一维热传导方程有限差分法的MATLAB实现》和曹刚教授的《一维偏微分方程的基本解》中，对偏微分方程的解得MATLAB实现问题进行过研究，但只停留在一维中，而实际中二维和三维的应用更加广泛。诸如粒子扩散或神经细胞的动作电位。也可以作为某些金融现象的模型，诸如布莱克-斯科尔斯模型与Ornstein-uhlenbeck过程。热方程及其非线性的推广形式也被应用与影响分析。 在科学和技术发展过程中，科学的理论和科学的实验一直是两种重要的科学方法和手段。虽然这两种科学方法都有十分重要的作用，但是一些研究对象往往由于他们的特性（例如太大或太小，太快或太慢）不能精确的用理论描述或用实验手段来实现。自从计算机出现和发展以来，模拟那些不容易观察到的现象，得到实际应用所需要的数值结果，解释各种现象的规律和基本性质。 科学计算在各门自然科学和技术科学与工程科学中其越来越大的作用，在很多重要领域中成为不可缺少的重要工具。而科学与工程计算中最重要的内容就是求解科学研究和工程技术中出现的各种各样的偏微分方程或方程组。 解偏微分方程已经成为科学与工程计算的核心内容，包括一些大型的计算和很多已经成为常规的计算。为什么它在当代能发挥这样大的作用呢？第一是计算机本身有了很大的发展；第二是数值求解方程的计算法有了很大的发展，这两者对人们计算能力的发展都是十分重要的。 1.2问题现状 近三十年来，解偏微分方程的理论和方法有了很大的发展，而且在各个学科技术的领域中应用也愈来愈广泛，在我国，偏微分方程数值解法作为一门课程，不但在计算数学专业，而且也在其他理工科专业的研究生的大学生中开设。同时，求解热传导方程的数值算法也取得巨大进展，特别是有限差分法方面，此算法的特点是在内边界处设计不同于整体的格式，将全局的隐式计算化为局部的分段隐式计算。而且精度上更好。 目前，在欧美各国MATLAB的使用十分普及。在大学的数学、工程和科学系科，MATLAB

抛物形扩散方程的有限差分法及数值实例

偏微分方程数值解 所在学院：数学与统计学院 课题名称：抛物形扩散方程的有限差分法及数值实例学生姓名：向聘

抛物形扩散方程的有限差分法及数值实例 1.1抛物型扩散方程 抛物型偏微分方程是一类重要的偏微分方程。考虑一维热传导方程： 22(),0u u a f x t T t x ??=+<≤?? (1.1.1) 其中a 是常数，()f x 是给定的连续函数。按照初边值条件的不同给法，可将(1.1.1)的定解分为两类： 第一，初值问题（Cauchy 问题）：求足够光滑的函数()t x u ,，满足方程(1.1.1)和初始条件： ()()x x u ?=0,, ∞<<∞-x (1.1.2) 第二，初边值问题（也称混合问题）：求足够光滑的函数()t x u ,，满足方程(1.1.1)和初始条件： ()()x x u ?=0,, 0x l << (1.1.3) 及边值条件 ()()0,,0==t l u t u ， T t ≤≤0 (1.1.4) 假定()x f 和()x ?在相应的区域光滑，并且于()0,0，()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。 1.2抛物线扩散方程的求解 下面考虑如下热传导方程 22()(0.)(,)0(,0)()u u a f x t x u t u L t u x x ????=+????? ==??=??? (1.2.1) 其中，0x l <<，T t ≤≤0,a (常数)是扩散系数。 取N l h = 为空间步长，M T =τ为时间步长，其中N ，M 是自然数，用两族

平行直线jh x x j ==, ()N j ,,1,0 =和k t t k τ ==, ()M k ,,1,0 =将矩形域 G {}T t l x ≤≤≤≤=0;0分割成矩形网格。其中 (),j k x t 表示网格节点；h G 表示 网格内点（位于开矩形G 中的网格节点）的集合；h G 表示位于闭矩形G 中的网格节点的集合；h Γ表示h G -h G 网格边界点的集合。 k j u 表示定义在网点(),j k x t 处的待求近似解，N j ≤≤0，M k ≤≤0。 现在对方程进行差分近似： (一) 向前差分格式 =-+τ k j k j u u 111 2 2(())k k k j j j j j j u u u a f f f x h +--++= (1.2.2) ()j j j x u ??==0， k u 0=k N u =0 (1.2.3) 计算后得： 111(12)k k k k j j j j j u ru r u ru f τ++-=+-++ (1.2.4) 其中，2 a r h τ = ，1,,1,0-=N j ，1,,1,0-=M k 。 显然，这是一个四点显示格式，每一层各个节点上的值是通过一个方程组求解到的。方程组如下： 1000 121011000 232121000 3432310001121(12)(12)(12)(12)N N N N N u ru r u ru f u ru r u ru f u ru r u ru f u ru r u ru f ττττ----?=+-++?=+-++??=+-++? ???=+-++? (1.2.5) 若记 () T k N k k k u u u 1 21,,,-= u ，()()()()T N x x x 121,,,-=???? ，()()()()T N x f x f x f 121,,,-=τττ f 则显格式(1.2.4)可写成向量形式 10 ,0,1,,1 k k k M φ +?=+=-?=? u Au f u (1.2.6) 其中

【文献综述】热传导方程差分格式的收敛性和稳定性

文献综述 信息与计算科学 热传导方程差分格式的收敛性和稳定性在实际研究物理问题过程中, 往往能给出问题相应的数学表达式, 但是由于实际物理问题的复杂性, 它的解却一般不容易求出. 由此计算物理应运而生, 计算物理是以计算机为工具, 应用数学的方法解决物理问题的一门应用性学科, 是物理、数学和计算机三者结合的交叉性学科. 它产生于二战期间美国对核武器的研究, 伴随着计算机的发展而发展. 计算物理的目的不仅仅是计算, 而是要通过计算来解释和发现新的物理规律. 这一点它与传统的实验物理和理论物理并无差别, 所不同的只是使用的工具和方法. 计算物理早已与实验物理和理论物理形成三足鼎立之势, 甚至有人提出它将成为现代物理大厦的“栋梁”. 在一个物理问题中一个数值解往往比一个式子更直观, 更有价值. 在实际求解方程时, 除了一些特殊的情况下可以方便地求得其精确解外, 在一般情况下, 当方程或定解条件具有比较复杂的形式, 或求解区域具有比较复杂的形状时, 往往求不到, 或不易求到其精确解. 这就需要我们去寻找方程的近似解, 特别是数值近似解, 简称数值解. 这里主要研究的是热传导方程. 有限差分法是微分方程和积分微分方程数值解的方法. 其基本思想是把连续的定解区域用有限个离散点构成的网格来代替, 这些离散点称作网格的节点;把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似;把原方程和定解条件中的微商用差商来近似, 积分用积分和来近似, 于是原微分方程和定解条件就近似地代之以代数方程组, 即有限差分方程组, 解此方程组就可以得到原问题在离散点上的近似解. 然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解. 热传导的差分法是求解热传导方程的重要方法之一. 对于差分格式的的求解, 我们首先要关注差分格式的收敛性和稳定性. 对于一个微分方程建立的各种差分格式, 为了有实用意义, 一个基本要求是它们能够任意逼近微分方程, 即相容性要求. 一个差分格式是否有用, 就要看差分方程的精确解能否任意逼近微分方程的解, 即收敛性的概念. 此外, 还有一个重要的概念必须考虑, 即差分格式的稳定性. 因为差分格式的计

抛物型方程的计算方法

分类号：O241.82 本科生毕业论文（设计） 题目：一类抛物型方程的计算方法 作者单位数学与信息科学学院 作者姓名 专业班级2011级数学与应用数学创新2班 指导教师 论文完成时间二〇一五年四月

一类抛物型方程的数值计算方法 (数学与信息科学学院数学与应用数学专业2011级创新2班) 指导教师 摘要: 抛物型方程数值求解常用方法有差分方法、有限元方法等。差分方法是一种对方程直接进行离散化后得到的差分计算格式，有限元方法是基于抛物型方程的变分形式给出的数值计算格式.本文首先给出抛物型方程的差分计算方法，并分析了相应差分格式的收敛性、稳定性等基本理论问题.然后，给出抛物型方程的有限元计算方法及理论分析. 关键词：差分方法，有限元方法，收敛性，稳定性 Numerical computation methods for a parabolic equation Yan qian (Class 2, Grade 2011, College of Mathematics and Information Science) Advisor: Nie hua Abstract: The common methods to solve parabolic equations include differential method, finite element method etc. The main idea of differential method is to construct differential schemes by discretizing differential equations directly. Finite element scheme is based on the variational method of parabolic equations. In this article, we give some differential schemes for a parabolic equation and analyze their convergence and stability. Moreover, the finite element method and the corresponding theoretical analysis for parabolic equation are established. Key words: differential method, finite element method, convergence, stability

抛物形扩散方程的有限差分法及数值实例(完整资料).doc

此文档下载后即可编辑 偏微分方程数值解 所在学院：数学与统计学院 课题名称：抛物形扩散方程的有限差分法及数值实例 学生姓名：向聘 抛物形扩散方程的有限差分法及数值实例

1.1抛物型扩散方程 抛物型偏微分方程是一类重要的偏微分方程。考虑一维热传导方程： 22(),0u u a f x t T t x ??=+<≤?? (1.1.1) 其中a 是常数，()f x 是给定的连续函数。按照初边值条件的不同给法，可将(1.1.1)的定解分为两类： 第一，初值问题（Cauchy 问题）：求足够光滑的函数()t x u ,，满足方程(1.1.1)和初始条件： ()()x x u ?=0,, ∞<<∞-x (1.1.2) 第二，初边值问题（也称混合问题）：求足够光滑的函数()t x u ,，满足方程(1.1.1)和初始条件： ()()x x u ?=0,, 0x l << (1.1.3) 及边值条件 ()()0,,0==t l u t u ， T t ≤≤0 (1.1.4) 假定()x f 和()x ?在相应的区域光滑，并且于()0,0，()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。 1.2抛物线扩散方程的求解 下面考虑如下热传导方程

22()(0.)(,)0(,0)()u u a f x t x u t u L t u x x ????=+????? ==??=??? (1.2.1) 其中，0x l <<，T t ≤≤0,a (常数)是扩散系数。 取N l h = 为空间步长，M T = τ为时间步长，其中N ，M 是自然 数，用两族平行直线 jh x x j ==, ()N j ,,1,0Λ=和 k t t k τ ==, ()M k ,,1,0Λ=将矩形域G {}T t l x ≤≤≤≤=0;0分割成矩形网格。其中 (),j k x t 表示网格节点；h G 表示网格内点（位于开矩形G 中的网格 节点）的集合；h G 表示位于闭矩形G 中的网格节点的集合；h Γ表示h G -h G 网格边界点的集合。 k j u 表示定义在网点(),j k x t 处的待求近似解， N j ≤≤0，M k ≤≤0。 现在对方程进行差分近似： (一) 向前差分格式 =-+τ k j k j u u 111 2 2(()) k k k j j j j j j u u u a f f f x h +--++= (1.2.2) () j j j x u ??==0， k u 0 = k N u =0 (1.2.3) 计算后得： 111(12)k k k k j j j j j u ru r u ru f τ++-=+-++ (1.2.4) 其中，2 a r h τ = ，1,,1,0-=N j Λ，1,,1,0-=M k Λ。 显然，这是一个四点显示格式，每一层各个节点上的值是通过一个方程组求解到的。方程组如下：

一维导热方程 有限差分法 matlab实现

第五次作业（前三题写在作业纸上） 一、用有限差分方法求解一维非定常热传导方程，初始条件和边界条件见说明.pdf 文件，热扩散系数α=const ， 22T T t x α??=?? 1. 用Tylaor 展开法推导出FTCS 格式的差分方程 2. 讨论该方程的相容性和稳定性，并说明稳定性要求对求解差分方程的影响。 3. 说明该方程的类型和定解条件，如何在程序中实现这些定解条件。 4. 编写M 文件求解上述方程，并用适当的文字对程序做出说明。（部分由网络搜索得到，添加，修改后得到。） function rechuandaopde %以下所用数据，除了t 的范围我根据题目要求取到了20000，其余均从pdf 中得来 a=0.00001;%a 的取值 xspan=[0 1];%x 的取值范围 tspan=[0 20000];%t 的取值范围 ngrid=[100 10];%分割的份数，前面的是t 轴的，后面的是x 轴的 f=@(x)0;%初值 g1=@(t)100;%边界条件一 g2=@(t)100;%边界条件二 [T,x,t]=pdesolution(a,f,g1,g2,xspan,tspan,ngrid);%计算所调用的函数 [x,t]=meshgrid(x,t); mesh(x,t,T);%画图，并且把坐标轴名称改为x ，t ，T xlabel('x') ylabel('t') zlabel('T') T%输出温度矩阵 dt=tspan(2)/ngrid(1);%t 步长 h3000=3000/dt;

h9000=9000/dt; h15000=15000/dt;%3000,9000,15000下，温度分别在T矩阵的哪些行T3000=T(h3000,:) T9000=T(h9000,:) T15000=T(h15000,:)%输出三个时间下的温度分布 %不再对三个时间下的温度-长度曲线画图，其图像就是三维图的截面 %稳定性讨论,傅里叶级数法 dx=xspan(2)/ngrid(2);%x步长 sta=4*a*dt/(dx^2)*(sin(pi/2))^2; if sta>0,sta<2 fprintf('\n%s\n','有稳定性') else fprintf('\n%s\n','没有稳定性') error end %真实值计算 [xe,te,Te]=truesolution(a,f,g1,g2,xspan,tspan,ngrid); [xe,te]=meshgrid(xe,te); mesh(xe,te,Te);%画图，并且把坐标轴名称改为xe，te，Te xlabel('xe') ylabel('te') zlabel('Te') Te%输出温度矩阵 %误差计算 jmax=1/dx+1;%网格点数 [rms]=wuchajisuan(T,Te,jmax) rms%输出误差

抛物型方程

前言 抛物型方程解的估计及其应用 1前言 数学物理方程主要指从物理学及其它各门自然科学、技术科学中所产生的偏微分方程（有时也包括积分方程、微分方程等），它们反映了有关的未知变量关于时间的导数和关于空间变量的导数之间的制约关系．连续介质力学、电磁学、量子力学等等方面的基本方程属于数学物理方程的范围．它以具有物理背景的偏微分方程（组）作为研究的主要对象．它与其他数学分支及物理、化学等自然科学和工程技术的很多领域都有着广泛的联系，因此，无论在历史上还是在今天的现实生活中，它对于推动数学理论的发展，加强理论与实际的联系，帮助人们认识世界和改造世界都起着重要大的作用． 微积分产生以后，人们就开始把力学中的一些问题，归结为偏微分方程进行研究．早在18世纪初，人们已经将弦线振动问题归结为弦振动方程，并探讨了它的解法．随后，人们又陆续了解了流体的运动、弹性体的平衡和振动、热传导、电磁相互作用、原子核和电子的相互作用、化学反应过程等等自然现象的基本规律，把它们写成偏微分方程的形式，并且求出了典型问题的解答，从而能通过实践，验证这些基本规律的正确性，显示了数学物理方程对于认识自然界基本规律的重要性．有了基本规律，人们还要利用这些基本规律来研究复杂的自然现象和解决复杂的工程技术问题，这就需要求出数学物理方程中许多特定问题的解答，随着计算机的出现及计算技术的发展，即使是相当复杂的问题，也可以计算出足够精确的数值来，这对于预测自然现象的变化(如气象预报)和进行各种工程设计（如机械强度的计算）都有着很重要的作用． 在研究数学物理方程的同时，人们对偏微分方程的性质也了解得越来越多、越来越深入，形成数学中的一门重要分支——偏微分方程理论．它既有悠久的历史，又不断地更新着它的对象、内容和方法．它直接联系着众多自然现象和实际问题，不断地提出或产生需要解决的新课题和新方法．它所面临的数学问题多样而复杂，不断地促进着许多相关数学分支（如泛函分析、复变函数、微分几何、计算数学等）的发展，并从它们之间引进许多有力的解决问题的工具．因此，数学物理方程又是纯粹数学的

有限差分法求解抛物型方程说明

有限差分法求解抛物型方程 偏微分方程只是在一些特殊情况下，才能求得定解问题解的解析式，对比较复杂的问题要找到解的解析表达式是困难的，因此需采用数值方法来求解.有限差分法是一种发展较早且比较成熟的数值求解方法，只适用于几何形状规则的结构化网格.它在微分方程中用差商代替偏导数，得到相应的差分方程，通过解差分方程得到微分方程解的近似值.本章主要介绍有限差分法的基本思想，并给出一些具体的数值实例. §1 差分方法的基本思想 有限差分法把偏微分方程的求解区域划分为有限个网格节点组成的网格，主要采用Taylor 级数展开等方法，在每个网格节点上用有限差分近似公式代替方程中的导数，从而建立以网格节点上的函数值为未知数的代数方程组. 有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式.从差分的空间形式来考虑，可分为中心格式和逆风格式.考虑时间因子的影响，差分格式还可以分为显格式、隐格式和显隐交替格式等.目前常见的差分格式，主要是上述几种格式的组合，不同的组合构成不同的差分格式. 泰勒级数展开法对有限差分格式的分类和公式的建立起着十分重要的作用.下面采用泰勒展开式导出一个自变量系统的若干有限差分表达式. 首先考虑单变量函数()u x ，如图1把区域x 离散为一批结点，记 0()(), =0,1,2,i i u x u x ih u i =+= 图1 单变量函数离散化 函数()u x 在点i x 处的泰勒展开式为 23 ()()()()()2!3! i i i i i u x u x u x h u x u x h h h ''''''+=++ ++ （1） 或 23 ()()()()()2!3! i i i i i u x u x u x h u x u x h h h ''''''-=-+ -+ （2） 式(1)和(2)重新整理可得 2()()()()()2!3! i i i i i u x h u x u x u x u x h h h '''''+-'= --- （3）

抛物形扩散方程的有限差分法与数值实例

偏微分方程数值解 所在学院:数学与统计学院 课题名称:抛物形扩散方程的有限差分法及数值实例学生：向聘

抛物形扩散方程的有限差分法及数值实例 1.1抛物型扩散方程 抛物型偏微分方程是一类重要的偏微分方程。考虑一维热传导方程: 2 u a 2 x 其中a 是常数，f(x)是给定的连续函数。按照初边值条件的不同给法，可将(1.1.1) 的定解分为两类： 第一，初值问题(Cauchy 问题)：求足够光滑的函数u x, t ，满足方程(1.1.1) 和初始条件： u x, 0 x , x (1.1.2) 第二，初边值问题(也称混合问题)：求足够光滑的函数u x,t ，满足方程(1.1.1) 和初始条件： u x, 0 x , 0 x l (1.1.3) 及边值条件 u0,t ul,t 0， 0 t T (1.1.4) 假定f X 和 x 在相应的区域光滑，并且于0,0，l,0两点满足相容条件， 则上述问题有唯一的充分光滑的解。 1.2抛物线扩散方程的求解 F 面考虑如下热传导方程 u(0.t) u(L,t) 0 u(x,0) (x) 其中，0 x l ，0 t T, a (常数)是扩散系数。 f(x),0 t T (1.1.1) (1.2.1)

取h秸为空间步长, M为时间步长，其中N，M是自然数，用两族

平行直线x X j jh , j 0,1, , N 和t t k k , k 0,1, , M 将矩形域 G 0 x l; 0 t T 分割成矩形网格。其中 X j ,t k 表示网格节点；G h 表示 网格点（位于开矩形G 中的网格节点）的集合；G h 表示位于闭矩形G 中的网格 节点的集合；h 表示G h - G h 网格边界点的集合。 u k 表示定义在网点x j ,t k 处的待求近似解，0 现在对方程进行差分近似: （一） 向前差分格式 若记 则显格式（1.2.4）可写成向量形式 U 其中 U j j X j , k k U 0 =U N =0 (1.2.3) 计算后得： 1 U j rU j 1 (1 2r)U k k r V 1 f j (1.2.4) 其中，r 各，j 0,1, , N 1，k h 0,1, ,M 1。 显然, 这是一个四点显示格式， 每一层各个节点上的值是通过一个方程组求 k 1 k U j U j h 2 k k k U j 1 2U j U j 1 a (f j f (X j )) 解到的。方程组如下: (1.2.2) 0 rU 2 (1 2r )U 0 0 叫 f 1 0 rU a (1 2r)U 0 rU 1 f 2 (1 0 rU 4 2r)U a rU 2 f 3 M 1 U N 1 0 「U N (1 2「)U N 1 「U N 2 f N (1.2.5) k k k k T u U 1 ,U 2, , U N 1 X 1 , X 2 , T X N 1 ， f f X ! , f X 2 , T X N 1 1 u Au 0,1,L ,M 1 (1.2.6)

有限差分法解偏微分方程

有限差分法解偏微分方程综述 绪论 有限元方法最早应用于结构力学，后来随着计算机的发展慢慢用于流体力学的数值模拟。在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的近似解构成。 有限差分方法(FDM)是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Taylor 级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。 对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。 从差分的空间形式来考虑，可分为中心格式和逆风格式。 考虑时间因子的影响，差分格式还可以分为显格式、隐格式、显隐交替格式等。 目前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。差分方法主要适用于有结构网格，网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。构造差分的方法有多种形式， 目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。 有限元方法的基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。采用不同的权函数和插值函数形式，便构成不同的有限元方法。 有限差分法求解偏微分方程 在采用数值计算方法求解偏微分方程时，若将每一处导数由有限差分近似公式替代，从而把求解偏微分方程的问题转换成求解代数方程的问题，即所谓的有限差分法。有限差分法求解偏微分方程的步骤如下： 1、区域离散化，即把所给偏微分方程的求解区域细分成由有限个格点组成的网格； 2、近似替代，即采用有限差分公式替代每一个格点的导数； 3、逼近求解。换而言之，这一过程可以看作是用一个插值多项式及其微分来代替偏微分方程的解的过程 有限差分法的应用 抛物型方程的差分方法 1. 简单差分法

抛物型方程有限差分法

抛物型方程有限差分法 1. 简单差分法 考虑一维模型热传导方程 (1.1) )(22x f x u a t u +??=??，T t ≤<0 其中a 为常数。)(x f 是给定的连续函数。（1.1）的定解问题分两类： 第一，初值问题（Cauchy 问题）：求足够光滑的函数()t x u ,，满足方程（1.1）和初始条件： （1.2） ()()x x u ?=0,, ∞<<∞-x 第二，初边值问题（也称混合问题）：求足够光滑的函数()t x u ,，满足方程（1.1）和初始条件： ()13.1 ()()x x u ?=0,, l x l <<- 及边值条件 ()23.1 ()()0,,0==t l u t u ， T t ≤≤0 假定()x f 和()x ?在相应的区域光滑，并且于()0,0，()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。 现在考虑边值问题（1.1），（1.3）的差分逼近 取 N l h = 为空间步长，M T =τ为时间步长，其中N ，M 是自然数， jh x x j ==, ()N j ,,1,0 =； τk y y k ==, ()M k ,,1,0 = 将矩形域G {}T t l x ≤≤≤≤=0;0分割成矩形网格。其中 ()j i y x ,表示网格节点； h G 表示网格内点（位于开矩形G 中的网格节点）的集合； h G 表示位于闭矩形G 中的网格节点的集合； h Γ表示h G -h G 网格边界点的集合。 k j u 表示定义在网点()k i t x ,处的待求近似解，N j ≤≤0，M k ≤≤0。

注意到在节点()k i t x ,处的微商和差商之间的下列关系（(,)k j k j u u x t t t ???? ≡ ?????） ： 可得到以下几种最简差分格式 （一） 向前差分格式 ()24.1 ()j j j x u ??==0， k u 0=k N u =0 其中1,,1,0-=N j ，1,,1,0-=M k 。取2h a r τ = 为网比，则进一步有 ()14.1' 1+k j u =k j ru 1++()r 21-k j u +k j ru 1-+j f τ 此差分格式是按层计算：首先，令0=k ，得到 1j u =01+j ru +()r 21-0j u +0 1-j ru +j f τ 于是，利用初值() j j j x u ??==0和边值k u 0=k N u =0，可算出第一层的1 j u ， 1,,1,0-=N j 。再由()14.1'取1=k ，可利用1j u 和k u 0=k N u =0算出2j u ， 1,,1,0-=N j 。如此下去，即可逐层算出所有k j u （1,,1,0-=N j ， 1,,1,0-=M k ） 。 由于第()1+k 层值可以通过第()k 层值直接得到，如此的格式称为显格式。并 视k j u 为()k j t x u ,的近似值。 若记 () T k N k k k u u u 1 21,,,-= u ，()()()()T N x x x 121,,,-=???? ，()()()()T N x f x f x f 121,,,-=τττ f 则显格式()14.1' 可写成向量形式 其中 若记 那末截断误差 （1.5） ()=u R k j () ()[]k j k j h Lu t x u L -,1=()ττO t x t u r k j +??? ? ??????? ??--)~,~(2112122=()2 h O +τ。 其中(,)j k x t 是矩形11+-<

有限差分法求解偏微分方程

有限差分法求解偏微分方程 摘要：本文主要使用有限差分法求解计算力学中的系统数学模型，推导了有限差分法的理论基础，并在此基础上给出了部分有限差分法求解偏微分方程的算例验证了推导的正确性及操作可行性。 关键词：计算力学，偏微分方程，有限差分法 Abstract：This dissertation mainly focuses on solving the mathematic model of computation mechanics with finite-difference method. The theoretical basis of finite-difference is derived in the second part of the dissertation, and then I use MATLAB to program the algorithms to solve some partial differential equations to confirm the correctness of the derivation and the feasibility of the method. Key words：Computation Mechanics, Partial Differential Equations, Finite-Difference Method

1 引言 机械系统设计常常需要从力学观点进行结构设计以及结构分析，而这些分析的前提就是建立工程问题的数学模型。通过对机械系统应用自然的基本定律和原理得到带有相关边界条件和初始条件的微分积分方程，这些微分积分方程构成了系统的数学模型。 求解这些数学模型的方法大致分为解析法和数值法两种，而解析法的局限性众所周知，当系统的边界条件和受载情况复杂一点，往往求不出问题的解析解或近似解。另一方面，计算机技术的发展使得计算更精确、更迅速。因此，对于绝大多数工程问题，研究其数值解法更具有实用价值。对于微分方程而言，主要分为差分法和积分法两种，本论文主要讨论差分法。 2 有限差分法理论基础 2.1 有限差分法的基本思想 当系统的数学模型建立后，我们面对的主要问题就是微分积分方程的求解。基本思想是用离散的只含有限个未知量的差分方程组去近似地代替连续变量的微分方程和定解条件，并把差分方程组的解作为微分方程定解问题的近似解。将原方程及边界条件中的微分用差分来近似，对于方程中的积分用求和或及机械求积公式来近似代替，从而把原微分积分方程和边界条件转化成差分方程组。有限差分法求解偏微分方程的步骤主要有以下几步： ?区域离散，即把所给偏微分方程的求解区域细分成由有限个格点组成的网格， 这些离散点称作网格的节点； ?近似替代，即采用有限差分公式替代每一个格点的导数； ?逼近求解，换而言之，这一过程可以看作是用一个插值多项式及其微分来代 替偏微分方程的解的过程。 从原则上说，这种方法仍然可以达到任意满意的计算精度。因为方程的连续数值解可以通过减小独立变量离散取值的间格，或者通过离散点上的函数值进行插值计算来近似得到。理论上，当网格步长趋近于零时，差分方程组的解应该收敛于精确解，但由于机器字节的限制，网格步长不可能也没有必要取得无限小，

热传导方程及其定解问题的导出

第一章 热传导方程 本章介绍最典型的抛物型方程—热传导方程,在研究热传导,扩散等物理现象时都会 遇到这类方程. §1 热传导方程及其定解问题的导出 热传导方程的导出 物理模型 在三维空间中,考虑一均匀,各向同性的物体Ω,假定它内部有热源,并且与周围介质有热交换,需要来研究物体内部温度的分布和变化. 以函数),,,(t z y x u 表示物体Ω在位置),,(z y x 及时刻t 的温度.物体内部由于各部分温度不同,产生热量的传递,它们遵循能量守恒定律. 能量守恒定律 物体内部的热量的增加等于通过物体的边界流入的热量与由物体内部的热源所生成的热量的总和. 在物体Ω内任意截取一块D .现在时段],[21t t 上对D 使用能量守恒定律. 设),,,(t z y x u u =是温度(度),c 是比热(焦耳∕度·千克),ρ是密度(千克/米3 ), q ρ 是 热流密度(焦耳/秒·米2 ),0f 是热源强度(焦耳/千克·秒). 注意到在dt 时段内通过D 的边界D ?上小块dS 进入区域D 的热量为dSdt n q ρρ?-(n ρ 是 D ?的外法向),从而由能量守恒律,我们有 ,)||(21 21 120??????????+?-=-?==t t D t t D D t t t t dxdydz f dt ds n q dt dxdydz u u c ρρρ ρ （） 大家知道,热量流动的原因是因为在物体内部存在温差.依据传热学中的傅立叶实验定律,在一定条件下,热流向量与温度梯度成正比 ,u k q ?-=ρ （梯度? ?? ? ????????==?z u y u x u gradu u ,,） （） 这里负号表明热量是由高温向低温流动,k 是物体的导热系数. ,n u k n u k n q ρρρρ??-=??-=? 从而式可改写为

二维抛物方程的有限差分法

二维抛物方程的有限差分法 摘要 二维抛物方程是一类有广泛应用的偏微分方程，由于大部分抛物方程都难以求得解析解，故考虑采用数值方法求解。有限差分法是最简单又极为重要的解微分方程的数值方法。本文介绍了二维抛物方程的有限差分法。 首先，简单介绍了抛物方程的应用背景，解抛物方程的常见数值方法，有限差分法的产生背景和发展应用。讨论了抛物方程的有限差分法建立的基础，并介绍了有限差分方法的收敛性和稳定性。其次，介绍了几种常用的差分格式，有古典显式格式、古典隐式格式、Crank-Nicolson隐式格式、Douglas差分格式、加权六点隐式格式、交替方向隐式格式等，重点介绍了古典显式格式和交替方向隐式格式。进行了格式的推导，分析了格式的收敛性、稳定性。并以热传导方程为数值算例，运用差分方法求解。通过数值算例，得出古典显式格式计算起来较简单，但稳定性条件较苛刻；而交替方向隐式格式无条件稳定。 关键词：二维抛物方程；有限差分法；古典显式格式；交替方向隐式格式

FINITE DIFFERENCE METHOD FOR TWO-DIMENSIONAL PARABOLIC EQUATION Abstract Two-dimensional parabolic equation is a widely used class of partial differential equations. Because this kind of equation is so complex, we consider numerical methods instead of obtaining analytical solutions. finite difference method is the most simple and extremely important numerical methods for differential equations. The paper introduces the finite difference method for two-dimensional parabolic equation. Firstly, this paper introduces the background and common numerical methods for Parabolic Equation, Background and development of applications. Discusses the basement for the establishment of the finite difference method for parabolic equation And describes the convergence and stability for finite difference method.Secondly, Introduces some of the more common simple differential format,for example, the classical explicit scheme, the classical implicit scheme, Crank-Nicolson implicit scheme, Douglas difference scheme, weighted six implicit scheme and the alternating direction implicit format. The paper focuses on the classical explicit scheme and the alternating direction implicit format. The paper takes discusses the derivation convergence,and stability of the format . The paper takes And the heat conduction equation for the numerical example, using the differential method to solve. Through numerical examples, the classical explicit scheme is relatively simple for calculation, with more stringent stability conditions; and alternating direction implicit scheme is unconditionally stable. Keywords:Two-dimensional Parabolic Equation; Finite-Difference Method; Eclassical Explicit Scheme; Alternating Direction Implicit Scheme