Chapter 3 有限差分法
3.1 泰勒展开式(Taylor Series)
在数学中,泰勒级数(英語:Taylor series)用无限项连加式——级数来表示一个函数,这些相加的项由函数在某一点的导数求得。泰勒级数是以于1715年发表了泰勒公式的英國数学家布魯克·泰勒(Sir Brook Taylor)来命名的。
泰勒展开式的基本形式:
f(x)=n∑k=0f(n)(0)n!(x)n
根据泰勒展开式,通过f(x)和其任意阶的导数,可以获得任意Δx值下的函数值f(x+Δx) ,即:
f(x+Δx)=f(x)+f′(x)1!Δx+f″(x)2!Δx2+f‴(x)3!Δx3+⋯+f(n)(x)n!(Δx)n
或者在−Δx位置,可写为:
f(x−Δx)=f(x)−f′(x)1!Δx+f″(x)2!Δx2−f‴(x)3!Δx3+⋯+f(n)(x)n!(−Δx)n
以上公式也可以写为: ui+1=ui+u′i1!Δx+u″i2!Δx2+u‴i3!Δx3+⋯+u(n)in!(Δx)n
ui−1=ui−u′i1!Δx+u″i2!Δx2−u‴i3!Δx3+⋯+u(n)in!(−Δx)n
在此,我们需要引入截断误差(truncation error)的概念,数学表达为O()。 O(2)和O(3)分别表示为在泰勒展开式上的二阶和三阶导数上的误差。截取误差的阶数越高,误差越小。
O(1)=u″i2!Δx2+u‴i3!Δx3+u(4)i4!Δx4+⋯+u(n)in!(Δx)n
O(2)=u‴i3!Δx3+⋯+u(n)in!(Δx)n O(3)=u(4)i4!Δx4+⋯+u(n)in!(Δx)n
何种情况下,泰勒展开式的截断误差为0?
- O(0)=0 时,意味着:f(x+Δx)=f(x)。则该函数为f(x)=C, C为常数。 如图:
- O(1)=0 时,意味着:f(x+Δx)=f(x)+f′(x)⋅Δx。则该函数为f(x)=ax+b。 如图:
如何依据泰勒级数,得到函数的一阶和二阶导数?
3.1.1 一阶导数
3.1.1.1 向前估计 (Forward Approximation)
采纳一阶截断误差,我们可将公式(3.2)写为: f(x+Δx)=f(x)+f′(x)1!Δx+O(1) 则: f′(x)=f(x+Δx)−f(x)Δx 或者 u′i=ui+1−uiΔx
注:公式(3.6)隐含了O(1)的误差。
3.1.1.2 向后估计 (Backward Approximation)
采纳一阶截断误差,我们可将公式(3.3)写为: f(x−Δx)=f(x)−f′(x)1!Δx+O(1) 则: f′(x)=f(x)−f(x−Δx)Δx 或者 u′i=ui−ui−1Δx
注:公式(3.7)隐含了O(1)的误差。
3.1.1.3 中心估计 (Central Approximation)
中心估计算法中,我们将从公式(3.2)减去公式(3.3),得到: f(x+Δx)−f(x−Δx)=0+2∗f′(x)1!Δx+0+2∗f‴(x)3!+...
截断误差由以上公式右边的第四项(三阶导数)开始,则该公式的截取误差为O(2),即二阶精度的截取误差,公式表达为: f(x+Δx)−f(x−Δx)=0+2∗f′(x)1!Δx+0+O(2) 可得到二阶精度的一阶导数的中心估计: f′(x)=f(x+Δx)−f(x−Δx)2Δx 或者 u′i=ui+1−ui−12Δx
公式(3.8)隐含了O(2)的误差,同时(3.6)和(3.7)都隐含了O(1)的误差。
3.1.2 二阶导数
f(x+Δx)+f(x−Δx)=2⋅f(x)+0+2⋅f″(x)2!Δx2+0+2⋅f(4)(x)4!Δx4+⋯
公式(3.9)来自公式 (3.2)和(3.3)的相加,三阶导数项在相加过程中为零,因此我们截取其三阶截取误差,则公式(3.9)可写为: f(x+Δx)+f(x−Δx)=2f(x)+f″(x)Δx2+O(3)
根据公式(3.10),我们可获得函数f(x)在x位置的二阶导数为:
f″(x)=f(x+Δx)−2f(x)+f(x−Δx)Δx2+O(3)
当移除其三阶截断误差O(3)后,我们得到近似的二阶导数:
f″(x)≈1Δx(f(x+Δx)−f(x)Δx+f(x)−f(x−Δx)Δx)≈f(x+Δx)−2f(x)+f(x−Δx)Δx2
将公式一般化,我们可写为以下形式:
u″i≈1Δx(ui+1−uiΔx−ui−ui−1Δx)≈ui+1−2ui+ui−1Δx2
3.2 构建数值方法
Example 3.1 一根100cm长的铁棍,初始温度25 ∘C,在其左右两边分别持续施加100∘C和50∘C的温度。 求解:任意时刻铁棍的温度分布。
参考信息:
问题描述
空间微分,如图。
物理定理
能量守恒:
能量变化 = 能量流入 - 能量流出 ΔE=Qin−Qout控制方程
通量=量单位时间⋅单位面积
- k - 热传导率[Wm−1K−1]。
- c - 比热容(specific heat capacity)[Jkg−1K−1]。
- ρ - 密度[kgm−3]。
- A - 截面积[m2]。
- D - 热力学扩散度(Thermal diffusivity)D=kρc [m2s−1]。
ρ∗c∗Δx∗A∗Δu=qin∗A∗Δt−qout∗A∗Δt 两边同时除以ρcΔxA,得到 ΔuΔt=1ρcqin−qoutΔx 以上公式当Δx趋近与0,Δt趋近于0时,得到微分形式: ∂u∂t=1ρc∂q∂x
q=k∂u∂x 则得到其控制方程: ∂u∂t=1ρc∂q∂x=1ρck∂u∂x∂x=D∂∂x(∂u∂x)
令D=kρc,单位[$ m s^{-2}$],则最终控制方程写为
∂u∂t=D∂2u∂x2
时空离散化
由一阶泰勒级数可知,控制方程(3.13)左边可写为: ∂u∂t=Duti−ut−1iΔt+O(1) 控制方程(3.13)左边可写为: ∂2u∂x2=ut−1i+1−2ut−1i+ut−1i−1Δx+O(3) 或者 ∂2u∂x2=uti+1−2uti+uti−1Δx+O(3)
此时,方程左边在时间尺度上具有一阶截取误差O(1),方程右边在空间尺度上具有三阶截断误差O(3)。省去误差项,离散化后控制方程写为: uti−ut−1iΔt=Dut−1i+1−2ut−1i+ut−1i−1Δx2 或者 uti−ut−1iΔt=Duti+1−2uti+uti−1Δx2
求解
3.3 显式求解法
显式求解法以公式(3.14)作为起点,该公式可变形为: uti−ut−1i=DΔtΔx2(ut−1i+1−2ut−1i+ut−1i−1) 令α=DΔtΔx2, β=1−2α,整理以上公式可得:
uti−ut−1i=α(ut−1i+1−2ut−1i+ut−1i−1)uti=αut−1i+1+(1−2α)ut−1i+αut−1i−1uti=αut−1i+1+βut−1i+αut−1i−1
将以上公式应用于离散点上,
点号 i | 公式 |
---|---|
1 | 边界条件: ut1=U0 |
2 | ut2=αut−13+βut−12+αut−11 |
3 | ut3=αut−14+βut−13+αut−12 |
4 | ut4=αut−15+βut−14+αut−13 |
5 | ut5=αut−16+βut−15+αut−14 |
… | … |
i-1 | … |
i | uti=αut−1i+1+βut−1i+αut−1i−1 |
i+1 | … |
… | … |
n-2 | utn−2=αut−1n−1+βut−1n−2+αut−1n−3 |
n-1 | utn−1=αut−1n+βut−1n−1+αut−1n−2 |
n | 边界条件: utn=UL |
由此我们得到n个算式,可转换为矩阵形式:
[ut1ut2ut3ut4…uti…utn−3utn−2utn−1utn]=[1000…0000α β α0…00000α β α…0000………………………00…α β α…00………………………0000…α β α00000…0α β α0000…0001][U0ut−12ut−13ut−14…ut−1i…ut−1n−3ut−1n−2ut−1n−1UL]
更简洁的方式,可写为:
[u1u2u3u4…ui…un−3un−2un−1un]t=[1000…0000α β α0…00000α β α…0000………………………00…α β α…00………………………0000…α β α00000…0α β α0000…0001][U0u2u3u4…ui…un−3un−2un−1UL]t−1
下一时刻(t)的变量组成的向量x由一个矩阵[A]乘以已知的前一时刻(t−1)的向量b获得,即:。 x=[A]∗b 由已知变量的矩阵求解未知变量的方法,称为显式求解法。
3.4 隐式求解法
显式求解法以公式(3.15)作为起点,该公式可变形为: uti−ut−1i=DΔtΔx2(uti+1−2uti+uti−1) 令α=DΔtΔx2, β=−1−2α,整理以上公式可得:
uti−ut−1i=α(uti+1−2uti+uti−1)−ut−1i=αuti+1+(−1−2α)uti+αuti−1−ut−1i=αuti+1+βuti+αuti−1
将以上公式应用于离散点上,
点号 i | 公式 |
---|---|
1 | 边界条件: ut−11=U0 |
2 | −ut−12=αut3+βut2+αut1 |
3 | −ut−13=αut4+βut3+αut2 |
4 | −ut−14=αut5+βut4+αut3 |
5 | −ut−15=αut6+βut5+αut4 |
… | … |
i-1 | … |
i | −ut−1i=αuti+1+βuti+αuti−1 |
i+1 | … |
… | … |
n-2 | −ut−1n−2=αutn−1+βutn−2+αutn−3 |
n-1 | −ut−1n−1=αutn+βutn−1+αutn−2 |
n | 边界条件: ut−1n=UL |
由此我们得到n个算式,可转换为矩阵形式:
[1000…0000α β α0…00000α β α…0000………………………00…α β α…00………………………0000…α β α00000…0α β α0000…0001][U0ut2ut3ut4…uti…utn−3utn−2utn−1UL]=−[ut−11ut−12ut−13ut−14…ut−1i…ut−1n−3ut−1n−2ut−1n−1ut−1n]
更简洁的方式,可写为:
[1000…0000α β α0…00000α β α…0000………………………00…α β α…00………………………0000…α β α00000…0α β α0000…0001][U0u2u3u4…ui…un−3un−2un−1UL]t=−[u1u2u3u4…ui…un−3un−2un−1un]t−1
一个矩阵[A]乘以下一时刻(t)的变量组成的向量x等于已知的前一时刻(t−1)的向量b,求解该方程则可得到x的值,数学表达为:
[A]∗x=b
通用的求解法为 x=[A]−1∗b
通过已知变量、未知变量和矩阵组成的公式或函数来求解未知变量的过程,称为隐式求解法。
3.5 编程求解
显式求解法
#' Problem: 1D Heat Transfer
#' governing Eqn: du/dt = k/r/c * (dd u / d x^2)
#' wiki: https://en.wikipedia.org/wiki/Thermal_diffusivity
#' wiki: https://en.wikipedia.org/wiki/Heat_equation
#' BC U0 = 100, UL=50
#' IC uic = 25
#' X = c(0, 1)
#' D = 23 mm2/s = 2.3e-5 m^2/s
#' DX = 0.01 m
#' DT = 10 s
#' Time = 0 to 1e6 s
HT.explicit <- function(
U0=100, UL=50, uic = 25,
X = 1,
DX= 0.1,
DT = 1,
DD = 2.3e-5,
Tmax = 1e5,
epsilon = 1e-4, bc2 = NULL, ignore.cfl = FALSE, plot = TRUE
){
T0 = 0
tt = seq(T0, Tmax, DT)
NT = length(tt)
xx = seq(0, X, DX)
NX = X / DX + 1
alpha = DD * DT / (DX * DX)
beta = 1-2*alpha
CFL = DD * DT / (DX * DX)
print(CFL)
if(!ignore.cfl){
if(CFL >=1 ){
stop()
}
}
#' =========================================
x0 = rep(uic, NX)
#' =========================================
ylim = c(min(x0), max(U0, UL))
xlim=c(0, X)
if(plot){
plot(xx, x0, type='b', col=2, lwd=3, ylim=ylim, xlim=xlim, xlab=xlab, ylab=ylab)
grid()
lines(x=c(1,1) * 0, y=c(min(x0), U0), lwd=2, col=3, type='b')
lines(x=c(1,1) * X, y=c(min(x0), UL), lwd=2, col=3, type='b')
text(x=X/2 , y = uic + diff(ylim)*.051, 'Initial condition', font=2)
text(x=X * 0.05 , y = U0, 'BC 1', font=2)
text(x=X * 0.95 , y = UL, 'BC 2', font=2)
}
#' =========================================
mat = matrix(0, nrow = NX, ncol = NX)
for(i in 1:NX){
for(j in 1:NX){
if(i==j){
mat[i, j] = beta
}
if(i+1 == j | i-1 == j ){
mat[i, j] = alpha
}
}
}
mat[1,]=c(1, rep(0, NX-1))
mat[NX,]=c(rep(0, NX-1), 1)
xm = matrix(NA, nrow=NX, ncol=NT)
vs = cbind(rep(0, NX))
# vs[NX/2] = ss
b=bx = cbind(x0)
for(i in 1:NT){
if(i>1){
bx = mat %*% b + vs * DT
if(any(is.nan(bx))) { break }
if(mean(abs(b-bx)) < epsilon) { break }
}
bx[1] = U0
bx[NX] = UL
xm[, i]=bx
b = bx
}
NT = i
xm=xm[, 1:NT]
# message('CFL value = ', CFL)
# message('Total Timesteps (dt * nt)= ', DT, ' * ', NT)
yy = xm; yy[abs(yy)>1e20] = NA
ret = list('x' = xx,
'time' = tt,
'u' = xm,
'CFL' = CFL,
'DT' = DT,
'NT' = NT,
'xlim' = xlim,
'ylim' = ylim)
return(ret)
}
plot1 <- function(x, nout = 20){
NT = x$NT
id=10^seq(0, log10(x$NT), length.out = nout)
col=colorspace::diverge_hcl(n=length(id)); lty=1
matplot(x=x$x, y=x$u[, id], type='l', ylim=x$ylim, xlim=x$xlim,
xlab=xlab, ylab=ylab, col=col, lty=lty)
legend('topright', paste0('T=', x$time[id]+1), col=col, lty=lty, bg='transparent')
mtext(text = paste('CFL =', x$CFL ), side=3, cex=1.5)
}
plot2 <- function(x){
NT = x$NT
id = c(2, 4, 6, 8); nid=length(id)
lty=1:nid; col=lty
matplot(t(x$u[id, ]), type='l', xlab=xlab, ylab=ylab, col=col, lty=lty); grid()
legend('bottomright',col=col, lty=lty, paste('Node', id))
}
隐式求解法
#' Problem: 1D Heat Transfer
#' governing Eqn: du/dt = k/r/c * (dd u / d x^2)
#' wiki: https://en.wikipedia.org/wiki/Thermal_diffusivity
#' wiki: https://en.wikipedia.org/wiki/Heat_equation
#' BC U0 = 100, UL=50
#' IC uic = 25
#' X = c(0, 1)
#' D = 23 mm2/s = 2.3e-5 m^2/s
#' DX = 0.01 m
#' DT = 10 s
#' Time = 0 to 1e6 s
HT.implicit <- function(
U0=100, UL=50, uic = 25,
X = 1,
DX= 0.1,
DT = 1,
DD = 2.3e-5,
Tmax = 1e5,
epsilon = 1e-4, bc2 = NULL, ignore.cfl = FALSE, plot = TRUE
){
T0 = 0
tt = seq(T0, Tmax, DT)
NT = length(tt)
xx = seq(0, X, DX)
NX = X / DX + 1
alpha = -DD * DT / (DX * DX)
beta = 1 + 2 * DD * DT / (DX * DX)
CFL = DD * DT / (DX * DX)
print(CFL)
if(!ignore.cfl){
if(CFL >=1 ){
stop()
}
}
#' =========================================
x0 = rep(uic, NX)
#' =========================================
ylim = c(min(x0), max(U0, UL))
xlim=c(0, X)
if(plot){
plot(xx, x0, type='b', col=2, lwd=3, ylim=ylim, xlim=xlim, xlab=xlab, ylab=ylab)
grid()
lines(x=c(1,1) * 0, y=c(min(x0), U0), lwd=2, col=3, type='b')
lines(x=c(1,1) * X, y=c(min(x0), UL), lwd=2, col=3, type='b')
text(x=X/2 , y = uic + diff(ylim)*.051, 'Initial condition', font=2)
text(x=X * 0.05 , y = U0, 'BC 1', font=2)
text(x=X * 0.95 , y = UL, 'BC 2', font=2)
}
#' =========================================
mat = matrix(0, nrow = NX, ncol = NX)
for(i in 1:NX){
for(j in 1:NX){
if(i==j){
mat[i, j] = beta
}
if(i+1 == j | i-1 == j ){
mat[i, j] = alpha
}
}
}
mat[1,]=c(1, rep(0, NX-1))
mat[NX,]=c(rep(0, NX-1), 1)
xm = matrix(NA, nrow=NX, ncol=NT)
vs = cbind(rep(0, NX))
# vs[NX/2] = ss
b=bx=cbind(x0)
for(i in 1:NT){
if(i>1){
bx = solve(mat, b)
if(any(is.nan(bx))) { break }
if(mean(abs(b-bx)) < epsilon) { break }
}
bx[1] = U0
bx[NX] = UL
xm[, i]=bx
b = bx
}
NT = i
xm=xm[, 1:NT]
# message('CFL value = ', CFL)
# message('Total Timesteps (dt * nt)= ', DT, ' * ', NT)
yy = xm; yy[abs(yy)>1e20] = NA
ret = list('x' = xx,
'time' = tt,
'u' = xm,
'CFL' = CFL,
'DT' = DT,
'NT' = NT,
'xlim' = xlim,
'ylim' = ylim)
return(ret)
}
plot1 <- function(x, nout = 20){
NT = x$NT
id=10^seq(0, log10(x$NT), length.out = nout)
col=colorspace::diverge_hcl(n=length(id)); lty=1
matplot(x=x$x, y=x$u[, id], type='l', ylim=x$ylim, xlim=x$xlim,
xlab=xlab, ylab=ylab, col=col, lty=lty)
legend('topright', paste0('T=', x$time[id]+1), col=col, lty=lty, bg='transparent')
mtext(text = paste('CFL =', x$CFL ), side=3, cex=1.5)
}
plot2 <- function(x){
NT = x$NT
id = c(2, 4, 6, 8); nid=length(id)
lty=1:nid; col=lty
matplot(t(x$u[id, ]), type='l', xlab=xlab, ylab=ylab, col=col, lty=lty); grid()
legend('bottomright',col=col, lty=lty, paste('Node', id))
}
3.6 显式与隐式求解法对比
3.6.1 CFL条件
source("Code/ch03/ch3_HeatTransferIm.R")
xlab ='Distance (m)'
ylab = 'Temperature (C)'
x = HT.implicit(DX= 0.05, DT = 1000, U0=100, UL=50, uic = 25, ignore.cfl=TRUE, plot=FALSE)
## [1] 9.2
3.6.2 计算效率
source("Code/ch03/ch3_HeatTransferIm.R")
source("Code/ch03/ch3_HeatTransferEx.R")
Tmax = 1e4
t0 = Sys.time()
x=HT.explicit(DX= 0.025, DT = 5, U0=100, UL=50, uic = 25, X = 1, ignore.cfl=FALSE, plot=FALSE, Tmax=Tmax)
## [1] 0.184
t1 = Sys.time()
tu.ex = t1 - t0
t0 = Sys.time()
x=HT.implicit(DX= 0.025, DT = 5, U0=100, UL=50, uic = 25, X = 1, ignore.cfl=FALSE, plot=FALSE, Tmax=Tmax)
## [1] 0.184
## Time for implicit =0.0979740619659424
## Time for explicit =0.06929612159729
3.7 二维有限差分
Example 3.2 sdhdt=kxB∗d2hdx2+kyB∗d2hdy2+ss
令Dx=kxBs和Dy=kyBs。
公式推导:
右边: ∂u∂t=ut+1i,j−uti,jΔt
左边:
Dx∂2u∂x2+Dy∂2u∂y2=Dxuti+1,j−2uti,j+uti−1,jΔx2+Dyuti,j+1−2uti,j+uti,j−1Δy2
控制方程离散化后得到: ut+1i,j−uti,jΔt=Dxuti+1,j−2uti,j+uti−1,jΔx2+Dyuti,j+1−2uti,j+uti,j−1Δy2
ut+1i,j−uti,j=DxΔtΔx2(uti+1,j−2uti,j+uti−1,j)+DxΔtΔy2(uti,j+1−2uti,j+uti,j−1)
令α=DxΔtΔx2, β=DxΔtΔy2, γ=1−2DxΔtΔx2−2DxΔtΔy2,公式变为:
ut+1i,j=αuti+1,j+βuti,j+1+γuti,j+βuti,j−1+αuti−1,j
假设x和y方向总长为Lx和Ly,沿两个方向的离散点数为Nx=Lx/Δx, Ny=Ly/Δy, N=Nx∗Ny。 矩阵形式可表达为: x=A∗b
x=[[u1,1…u1,Ny][u2,1…u2,Ny]…ui,j…[uNx,1…uNx,Ny]]t
A=[[100010001]0000[1000000β…αγα…β0000001]0000[1000000β…αγα…β0000001]0000[100010001]]
b=[[U1]Ny∗1[U2u2,1…U3]Ny∗1[U2ui,j…U3]Ny∗1[U4]Ny∗1]t−1
3.7.1 编程求解
显式求解
#' Problem: 1D Confined Aquifer
#' governing Eqn: du/dt = DDx * (dd u / d x^2) + DDy * (dd u / d y^2)
#'
diag.matrix <- function(id = c(-1, 0, 1), x = rep(1, length(id)), n = 3, def.val = 0){
val = matrix(x, ncol=length(id), nrow=1)
mat = matrix(def.val, n, n)
nid = length(id)
for(i in 1:n){
for(j in 1:n){
for(k in 1:nid){
if(i + id[k] == j){
mat[i,j] = val[k]
}
}
}
}
mat
}
toBC <- function(idl, x, val){
nbc = length(idl)
for(i in 1:nbc){
x[idl[[i]]] = val[i]
}
x
}
CA.Explicit <- function( bc1 = c(0, 0, 0,0),
bc2 = NULL,
uic = 25,
Lxy = c(1000, 1000),
Dxy = c(50,50),
DD = rep(23 ,2),
epsilon = 0.001,
DT = 25,
Tmax = 1e5,
ignore.cfl = FALSE, plot = TRUE){
DX=Dxy[1]; DY=Dxy[2];
tt = seq(0, Tmax, DT);
NT = length(tt)
xx = seq(0, Lxy[1], DX); NX = length(xx)
yy = seq(0, Lxy[2], DY); NY = length(yy)
# NX = Lxy[1] / DX + 1; NY = Lxy[2] / DY + 1
N = NX * NY
CFL.x = alpha = DD[1] * DT / (DX * DX)
CFL.y = beta = DD[2] * DT / (DY * DY)
gamma = 1 - 2 * alpha - 2 * beta
message('CFL value = (', CFL.x, '\t', CFL.y, ')')
if(!ignore.cfl){
if(CFL.x >=.5 | CFL.y >=.5){
stop()
}
}
mat = diag.matrix(id = c(-NY, -1, 0, 1, NY), n=N,
x=c(alpha, beta, gamma, beta, alpha), def.val = 0)
dmat = diag.matrix(id=0, x=1, n=N, def.val = 0)
idl = list(1:NY,
1+(1:NX - 1)*(NY),
(1:NX) * NY,
(NX-1)*(NY)+1:NY)
nbc = length(idl)
id.bc = sort(unique(unlist(idl)))
mat[id.bc, ] = dmat[id.bc,]
arr = array(NA, dim=c(NY,NX,NT))
# xm = matrix(NA, nrow=N, ncol=NT)
vs = cbind(rep(0, N))
# bc2=list(id=10, x=0.01)
# vs[bc2$id] = bc2$x
b=bx=cbind(rep(uic, N))
b = toBC(idl = idl, x=b, val=bc1)
for(i in 1:NT){
if(i>1){
bx = mat %*% b + vs * DT
if(any(is.nan(bx))) { break }
if(mean(abs(b-bx)) < epsilon) { break }
}
bx = toBC(idl = idl, x=bx, val=bc1)
arr[, , i] = matrix(bx, nrow = NY, ncol = NX)
b = bx
}
NT = i
arr=arr[,, 1:NT]
# message('Total Timesteps (dt * nt)= ', DT, ' * ', NT)
# yy = arr; yy[abs(yy)>1e20] = NA
ret = list('x' = xx,
'y' = yy,
'z' = arr,
'time' = tt,
'CFL' = c(CFL.x, CFL.y),
'DT' = DT,
'NT' = NT
)
return(ret)
}
plot.3d <- function(x, nr=3, nc=4, clim=NULL){
par(mfrow=c(nr, nc), mar=c(1,1,1,1))
idx = round(10^seq(0, log10(x$NT), length.out = nc*nr))
z=x$z;
z[ is.infinite(abs(z)) ] = NA
if(is.null(clim)){
clim = range(z, na.rm = TRUE)
}
for(i in idx ){
plot3D::persp3D(z=z[, , i], clim=clim,
colvar=z[, , i])
mtext(paste('T =', x$DT * i), side= 3, line=-1)
}
}
library(plot3D)
source("Code/ch03/ch3_ConfAq2DIm.R")
x = CA.Explicit(bc1 = c(10, 30,10, 10),
uic = 25,
Lxy = c(1000, 1000),
Dxy = c(50,50),
DD = rep(23 ,2),
epsilon = 0.001,
DT = 25,
Tmax = 1000
)
## CFL value = (0.23 0.23)
隐式求解
#' Problem: 1D Confined Aquifer
#' governing Eqn: du/dt = DDx * (dd u / d x^2) + DDy * (dd u / d y^2)
#'
diag.matrix <- function(id = c(-1, 0, 1), x = rep(1, length(id)), n = 3, def.val = 0){
val = matrix(x, ncol=length(id), nrow=1)
mat = matrix(def.val, n, n)
nid = length(id)
for(i in 1:n){
for(j in 1:n){
for(k in 1:nid){
if(i + id[k] == j){
mat[i,j] = val[k]
}
}
}
}
mat
}
toBC <- function(idl, x, val){
nbc = length(idl)
for(i in 1:nbc){
x[idl[[i]]] = val[i]
}
x
}
CA.Implicit <- function( bc1 = c(0, 0, 0,0),
bc2 = NULL,
uic = 25,
Lxy = c(1000, 1000),
Dxy = c(50,50),
DD = rep(23 ,2),
epsilon = 0.001,
DT = 25,
Tmax = 1e5,
ignore.cfl = FALSE, plot = TRUE){
DX=Dxy[1]; DY=Dxy[2];
tt = seq(0, Tmax, DT);
NT = length(tt)
xx = seq(0, Lxy[1], DX); NX = length(xx)
yy = seq(0, Lxy[2], DY); NY = length(yy)
# NX = Lxy[1] / DX + 1; NY = Lxy[2] / DY + 1
N = NX * NY
alpha = -DD[1] * DT / (DX * DX)
beta = -DD[2] * DT / (DY * DY)
CFL.x = DD[1] * DT / (DX * DX)
CFL.y = DD[2] * DT / (DY * DY)
gamma = 1 + 2 * DD[1] * DT / (DX * DX) + 2 * DD[2] * DT / (DY * DY)
message('CFL value = (', CFL.x, '\t', CFL.y, ')')
if(!ignore.cfl){
if(CFL.x >=.5 | CFL.y >=.5){
stop()
}
}
#' =========================================
x0 = rep(uic, NX)
mat = diag.matrix(id = c(-NY, -1, 0, 1, NY), n=N,
x=c(alpha, beta, gamma, beta, alpha), def.val = 0)
dmat = diag.matrix(id=0, x=1, n=N, def.val = 0)
idl = list(1:NY,
1+(1:NX - 1)*(NY),
(1:NX) * NY,
(NX-1)*(NY)+1:NY)
nbc = length(idl)
id.bc = sort(unique(unlist(idl)))
mat[id.bc, ] = dmat[id.bc,]
arr = array(NA, dim=c(NY,NX,NT))
vs = cbind(rep(0, N))
b=bx=cbind(rep(uic, N))
b = toBC(idl = idl, x=b, val=bc1)
for(i in 1:NT){
if(i>1){
bx = solve(mat, b + vs * DT)
if(any(is.nan(bx))) { break }
if(mean(abs(b-bx)) < epsilon) { break }
}
bx = toBC(idl = idl, x=bx, val=bc1)
arr[, , i] = matrix(bx, nrow = NY, ncol = NX)
b = bx
}
NT = i
arr=arr[,, 1:NT]
# message('Total Timesteps (dt * nt)= ', DT, ' * ', NT)
# yy = arr; yy[abs(yy)>1e20] = NA
ret = list('x' = xx,
'y' = yy,
'z' = arr,
'time' = tt,
'CFL' = c(CFL.x, CFL.y),
'DT' = DT,
'NT' = NT
)
return(ret)
}
plot.3d <- function(x, nr=3, nc=4, clim=NULL){
par(mfrow=c(nr, nc), mar=c(1,1,1,1))
idx = round(10^seq(0, log10(x$NT), length.out = nc*nr))
z=x$z;
z[ is.infinite(abs(z)) ] = NA
if(is.null(clim)){
clim = range(z, na.rm = TRUE)
}
for(i in idx ){
plot3D::persp3D(z=z[, , i], clim=clim,
colvar=z[, , i])
mtext(paste('T =', x$DT * i), side= 3, line=-1)
}
}
library(plot3D)
source("Code/ch03/ch3_ConfAq2DIm.R")
x = CA.Implicit(bc1 = c(10, 30,10, 10),
uic = 25,
Lxy = c(1000, 1000),
Dxy = c(50,50),
DD = rep(23 ,2),
epsilon = 0.001,
DT = 25,
Tmax = 1000
)
## CFL value = (0.23 0.23)
对比隐式与显式求解法的时间步长和效率:
rm(list=ls())
library(plot3D)
source("Code/ch03/ch3_ConfAq2DEx.R")
source("Code/ch03/ch3_ConfAq2DIm.R")
Tmax = 1e4
dx = dy = 50
t0 = Sys.time()
x1 = CA.Explicit(Tmax = Tmax, DT = 25,
Lxy = c(1000, 1000), Dxy = c(dx, dy),
ignore.cfl = TRUE)
## CFL value = (0.23 0.23)
t1 = Sys.time()
tu.ex = t1 - t0
t0 = Sys.time()
x2 = CA.Implicit(Tmax = Tmax, DT = 25,
Lxy = c(1000, 1000), Dxy = c(dx, dy),
ignore.cfl = TRUE)
## CFL value = (0.23 0.23)
t1 = Sys.time()
tu.im = t1 - t0
df = data.frame( 'Type' = c('Explicit', 'Implicit'),
'NT' = c(x1$NT, x2$NT),
'DT' = c(x1$DT, x2$DT),
'Time' = c(tu.ex, tu.im),
'CFL.X' = c(x1$CFL[1], x2$CFL[1]),
'CFL.Y' = c(x1$CFL[2], x2$CFL[2])
)
print(df)
## Type NT DT Time CFL.X CFL.Y
## 1 Explicit 401 25 3.523902 secs 0.23 0.23
## 2 Implicit 401 25 7.618441 secs 0.23 0.23