Chapter 3 有限差分法

3.1 泰勒展开式(Taylor Series)

在数学中,泰勒级数(英語:Taylor series)用无限项连加式——级数来表示一个函数,这些相加的项由函数在某一点的导数求得。泰勒级数是以于1715年发表了泰勒公式的英國数学家布魯克·泰勒(Sir Brook Taylor)来命名的。

泰勒展开式的基本形式:

f(x)=nk=0f(n)(0)n!(x)n

Taylor Series
Taylor Series

根据泰勒展开式,通过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!Δx2f(x)3!Δx3++f(n)(x)n!(Δx)n

以上公式也可以写为: ui+1=ui+ui1!Δx+ui2!Δx2+ui3!Δx3++u(n)in!(Δx)n

ui1=uiui1!Δx+ui2!Δx2ui3!Δx3++u(n)in!(Δx)n

在此,我们需要引入截断误差(truncation error)的概念,数学表达为O()O(2)O(3)分别表示为在泰勒展开式上的二阶和三阶导数上的误差。截取误差的阶数越高,误差越小

O(1)=ui2!Δx2+ui3!Δx3+u(4)i4!Δx4++u(n)in!(Δx)n

O(2)=ui3!Δ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 或者 ui=ui+1uiΔ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 或者 ui=uiui1Δx

注:公式(3.7)隐含了O(1)的误差。

向后估计
向后估计

3.1.1.3 中心估计 (Central Approximation)

中心估计算法中,我们将从公式(3.2)减去公式(3.3),得到: f(x+Δx)f(xΔx)=0+2f(x)1!Δx+0+2f(x)3!+...

截断误差由以上公式右边的第四项(三阶导数)开始,则该公式的截取误差为O(2),即二阶精度的截取误差,公式表达为: f(x+Δx)f(xΔx)=0+2f(x)1!Δx+0+O(2) 可得到二阶精度的一阶导数的中心估计: f(x)=f(x+Δx)f(xΔx)2Δx 或者 ui=ui+1ui12Δx

公式(3.8)隐含了O(2)的误差,同时(3.6)(3.7)都隐含了O(1)的误差。

中心估计
中心估计

3.1.2 二阶导数

我们将公式 (3.2)(3.3)相加,可得到:

f(x+Δx)+f(xΔx)=2f(x)+0+2f(x)2!Δx2+0+2f(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

将公式一般化,我们可写为以下形式:

ui1Δx(ui+1uiΔxuiui1Δx)ui+12ui+ui1Δx2

二阶导数估计
二阶导数估计

3.2 构建数值方法

Example 3.1 一根100cm长的铁棍,初始温度25 C,在其左右两边分别持续施加100C50C的温度。 求解:任意时刻铁棍的温度分布。

参考信息:

  1. 问题描述

    空间微分,如图。

  2. 物理定理

    能量守恒:
    能量变化 = 能量流入 - 能量流出 ΔE=QinQout

  3. 控制方程

    =

    • k - 热传导率[Wm1K1]。
    • c - 比热容(specific heat capacity)[Jkg1K1]。
    • ρ - 密度[kgm3]。
    • A - 截面积[m2]。
    • D - 热力学扩散度(Thermal diffusivity)D=kρc [m2s1]。

    ρcΔxAΔu=qinAΔtqoutAΔt 两边同时除以ρcΔxA,得到 ΔuΔt=1ρcqinqoutΔx 以上公式当Δx趋近与0,Δt趋近于0时,得到微分形式: ut=1ρcqx

    q=kux 则得到其控制方程: ut=1ρcqx=1ρckuxx=Dx(ux)

    D=kρc,单位[$ m s^{-2}$],则最终控制方程写为

    ut=D2ux2

  4. 时空离散化

    由一阶泰勒级数可知,控制方程(3.13)左边可写为: ut=Dutiut1iΔt+O(1) 控制方程(3.13)左边可写为: 2ux2=ut1i+12ut1i+ut1i1Δx+O(3) 或者 2ux2=uti+12uti+uti1Δx+O(3)

    此时,方程左边在时间尺度上具有一阶截取误差O(1),方程右边在空间尺度上具有三阶截断误差O(3)。省去误差项,离散化后控制方程写为: utiut1iΔt=Dut1i+12ut1i+ut1i1Δx2 或者 utiut1iΔt=Duti+12uti+uti1Δx2

  5. 求解

3.3 显式求解法

显式求解法以公式(3.14)作为起点,该公式可变形为: utiut1i=DΔtΔx2(ut1i+12ut1i+ut1i1)α=DΔtΔx2, β=12α,整理以上公式可得:

utiut1i=α(ut1i+12ut1i+ut1i1)uti=αut1i+1+(12α)ut1i+αut1i1uti=αut1i+1+βut1i+αut1i1

将以上公式应用于离散点上,

点号 i 公式
1 边界条件: ut1=U0
2 ut2=αut13+βut12+αut11
3 ut3=αut14+βut13+αut12
4 ut4=αut15+βut14+αut13
5 ut5=αut16+βut15+αut14
i-1
i uti=αut1i+1+βut1i+αut1i1
i+1
n-2 utn2=αut1n1+βut1n2+αut1n3
n-1 utn1=αut1n+βut1n1+αut1n2
n 边界条件: utn=UL

由此我们得到n个算式,可转换为矩阵形式:

[ut1ut2ut3ut4utiutn3utn2utn1utn]=[10000000α β α000000α β  α000000α β α000000α  β α000000α β α00000001][U0ut12ut13ut14ut1iut1n3ut1n2ut1n1UL]

更简洁的方式,可写为:

[u1u2u3u4uiun3un2un1un]t=[10000000α β α000000α β  α000000α β α000000α  β α000000α β α00000001][U0u2u3u4uiun3un2un1UL]t1

下一时刻(t)的变量组成的向量x由一个矩阵[A]乘以已知的前一时刻(t1)的向量b获得,即:。 x=[A]b 由已知变量的矩阵求解未知变量的方法,称为显式求解法

3.4 隐式求解法

显式求解法以公式(3.15)作为起点,该公式可变形为: utiut1i=DΔtΔx2(uti+12uti+uti1)α=DΔtΔx2, β=12α,整理以上公式可得:

utiut1i=α(uti+12uti+uti1)ut1i=αuti+1+(12α)uti+αuti1ut1i=αuti+1+βuti+αuti1

将以上公式应用于离散点上,

点号 i 公式
1 边界条件: ut11=U0
2 ut12=αut3+βut2+αut1
3 ut13=αut4+βut3+αut2
4 ut14=αut5+βut4+αut3
5 ut15=αut6+βut5+αut4
i-1
i ut1i=αuti+1+βuti+αuti1
i+1
n-2 ut1n2=αutn1+βutn2+αutn3
n-1 ut1n1=αutn+βutn1+αutn2
n 边界条件: ut1n=UL

由此我们得到n个算式,可转换为矩阵形式:

[10000000α β α000000α β  α000000α β α000000α  β α000000α β α00000001][U0ut2ut3ut4utiutn3utn2utn1UL]=[ut11ut12ut13ut14ut1iut1n3ut1n2ut1n1ut1n]

更简洁的方式,可写为:

[10000000α β α000000α β  α000000α β α000000α  β α000000α β α00000001][U0u2u3u4uiun3un2un1UL]t=[u1u2u3u4uiun3un2un1un]t1

一个矩阵[A]乘以下一时刻(t)的变量组成的向量x等于已知的前一时刻(t1)的向量b,求解该方程则可得到x的值,数学表达为:

[A]x=b

通用的求解法为 x=[A]1b

通过已知变量、未知变量和矩阵组成的公式或函数来求解未知变量的过程,称为隐式求解法

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
plot1(x, nout=10)

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
t1 = Sys.time()
tu.im = t1 - t0

message('Time for implicit =', tu.im)
## Time for implicit =0.0979740619659424
message('Time for explicit =', tu.ex)
## Time for explicit =0.06929612159729

3.7 二维有限差分

Example 3.2 sdhdt=kxBd2hdx2+kyBd2hdy2+ss

Dx=kxBsDy=kyBs

公式推导:

右边: ut=ut+1i,juti,jΔt

左边:

Dx2ux2+Dy2uy2=Dxuti+1,j2uti,j+uti1,jΔx2+Dyuti,j+12uti,j+uti,j1Δy2

控制方程离散化后得到: ut+1i,juti,jΔt=Dxuti+1,j2uti,j+uti1,jΔx2+Dyuti,j+12uti,j+uti,j1Δy2

ut+1i,juti,j=DxΔtΔx2(uti+1,j2uti,j+uti1,j)+DxΔtΔy2(uti,j+12uti,j+uti,j1)

α=DxΔtΔx2, β=DxΔtΔy2, γ=12DxΔtΔx22DxΔtΔy2,公式变为:

ut+1i,j=αuti+1,j+βuti,j+1+γuti,j+βuti,j1+αuti1,j

二维离散化格点
二维离散化格点

假设xy方向总长为LxLy,沿两个方向的离散点数为Nx=Lx/Δx, Ny=Ly/Δy, N=NxNy。 矩阵形式可表达为: x=Ab

x=[[u1,1u1,Ny][u2,1u2,Ny]ui,j[uNx,1uNx,Ny]]t

A=[[100010001]0000[1000000βαγαβ0000001]0000[1000000βαγαβ0000001]0000[100010001]]

b=[[U1]Ny1[U2u2,1U3]Ny1[U2ui,jU3]Ny1[U4]Ny1]t1

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)
plot.3d(x)

隐式求解

#' 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)
plot.3d(x)

对比隐式与显式求解法的时间步长和效率:

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