Merge branch 'advection_naaice' into 'main'
Merge heterogeneous diffusion See merge request naaice/tug!1
This commit is contained in:
Marco De Lucia
commit
79957da687
1
.gitignore
vendored
1
.gitignore
vendored
@ -8,3 +8,4 @@ compile_commands.json
|
||||
.cache/
|
||||
.ccls-cache/
|
||||
/iwyu/
|
||||
.Rhistory
|
||||
|
||||
@ -5,11 +5,12 @@
|
||||
#+OPTIONS: toc:nil
|
||||
|
||||
|
||||
* Diffusion in 1D
|
||||
* Homogeneous diffusion in 1D
|
||||
|
||||
** Finite differences with nodes as cells' centres
|
||||
|
||||
The 1D diffusion equation is:
|
||||
The 1D diffusion equation for spatially constant diffusion
|
||||
coefficients $\alpha$ is:
|
||||
|
||||
#+NAME: eqn:1
|
||||
\begin{align}
|
||||
@ -304,3 +305,137 @@ form:
|
||||
|
||||
|
||||
#+LATEX: \clearpage
|
||||
|
||||
* Heterogeneous diffusion
|
||||
|
||||
If the diffusion coefficient $\alpha$ is spatially variable, [[eqn:1]] can
|
||||
be rewritten:
|
||||
|
||||
#+NAME: eqn:hetdiff
|
||||
\begin{align}
|
||||
\frac{\partial C }{\partial t} & = \frac{\partial}{\partial x} \left(\alpha(x) \frac{\partial C }{\partial x} \right)
|
||||
\end{align}
|
||||
|
||||
** Discretization of the equation using chain rule
|
||||
|
||||
\noindent From the product rule for derivatives we obtain:
|
||||
|
||||
#+NAME: eqn:product
|
||||
\begin{equation}
|
||||
\frac{\partial}{\partial x} \left(\alpha(x) \frac{\partial C }{\partial x}\right) = \frac{\partial \alpha}{\partial x} \cdot \frac{\partial C}{\partial x} + \alpha \frac{\partial^2 C }{\partial x^2}
|
||||
\end{equation}
|
||||
|
||||
\noindent Using a spatially centred second order finite difference
|
||||
approximation at $x=x_i$ for both $\alpha$ and $C$, we have
|
||||
#+NAME: eqn:hetdiff_fd
|
||||
\begin{align}
|
||||
\frac{\partial \alpha}{\partial x} \cdot \frac{\partial C}{\partial x} +
|
||||
\alpha \frac{\partial^2 C }{\partial x^2} & \simeq \frac{\alpha_{i+1} - \alpha_{i-1}}{2\Delta x}\cdot\frac{C_{i+1} - C_{i-1}}{2\Delta x} + \alpha_i \frac{C_{i+1} - 2 C_{i} + C_{i-1}}{\Delta x^2} \\ \nonumber
|
||||
& = \frac{1}{\Delta x^2} \frac{\alpha_{i+1} - \alpha_{i-1}}{4}(C_{i+1} - C_{i-1}) + \frac{\alpha_{i}}{\Delta x^2}(C_{i+1}-2C_i+C_{i-1})\\ \nonumber
|
||||
& = \frac{1}{\Delta x^2} \left\{A C_{i+1} -2\alpha_i C_i + AC_{i-1})\right\}
|
||||
\end{align}
|
||||
\noindent having set
|
||||
\[ A = \frac{\alpha_{i+1}-\alpha_{i-1}}{4} + \alpha_i \]
|
||||
|
||||
\noindent In 2D the ADI scheme [[eqn:genADI]] with heterogeneous diffusion
|
||||
coefficients can thus be written:
|
||||
|
||||
#+NAME: eqn:genADI_het
|
||||
\begin{equation}
|
||||
\systeme{
|
||||
\displaystyle \frac{C^{t+1/2}_{i,j}-C^{t }_{i,j}}{\Delta t/2} = \displaystyle \frac{\partial}{\partial x} \left( \alpha^x_{i,j} \frac{\partial C^{t+1/2}_{i,j}}{\partial x}\right) + \frac{\partial}{\partial y} \left( \alpha^y_{i,j} \frac{\partial C^{t }_{i,j}}{\partial y}\right),
|
||||
\displaystyle \frac{C^{t+1 }_{i,j}-C^{t+1/2}_{i,j}}{\Delta t/2} = \displaystyle \frac{\partial}{\partial x} \left( \alpha^x_{i,j} \frac{\partial C^{t+1/2}_{i,j}}{\partial x}\right) + \frac{\partial}{\partial y} \left( \alpha^y_{i,j} \frac{\partial C^{t+1}_{i,j}}{\partial y}\right)
|
||||
}
|
||||
\end{equation}
|
||||
|
||||
\noindent We define for compactness $S_x=\frac{\Delta t}{2\Delta x^2}$
|
||||
and $S_y=\frac{\Delta t}{2\Delta y^2}$ and
|
||||
|
||||
#+NAME: eqn:het_AB
|
||||
\begin{equation}
|
||||
\systeme{
|
||||
\displaystyle A_{i,j} = \displaystyle \frac{\alpha^x_{i+1,j} -\alpha^x_{i-1,j }}{4} + \alpha^x_{i,j},
|
||||
\displaystyle B_{i,j} = \displaystyle \frac{\alpha^y_{i, j+1}-\alpha^y_{i ,j-1}}{4} + \alpha^y_{i,j}
|
||||
}
|
||||
\end{equation}
|
||||
|
||||
\noindent Plugging eq. ([[eqn:hetdiff_fd]]) into the first of equations
|
||||
([[eqn:genADI_het]]) - so called "sweep by x" - and putting all implicit
|
||||
terms (at time level $t+1/2$) on the left hand side we obtain:
|
||||
|
||||
#+NAME: eqn:sweepX_het
|
||||
\begin{equation}\displaystyle
|
||||
\begin{split}
|
||||
-S_x A_{i,j} C^{t+1/2}_{i+1,j} + (1 + 2S_x\alpha^x_{i,j})C^{t+1/2}_{i,j} - S_x A_{i,j}C^{t+1/2}_{i-1,j} = \\
|
||||
S_y B_{i,j} C^{t }_{i,j+1} + (1 - 2S_y\alpha^y_{i,j})C^{t }_{i,j} + S_y B_{i,j}C^{t }_{i,j-1}
|
||||
\end{split}
|
||||
\end{equation}
|
||||
|
||||
\noindent In the same way for the second of eq. [[eqn:genADI_het]] we have:
|
||||
|
||||
#+NAME: eqn:sweepY_het
|
||||
\begin{equation}\displaystyle
|
||||
\begin{split}
|
||||
-S_y B_{i,j} C^{t+1 }_{i,j+1} + (1 + 2S_y\alpha^y_{i,j})C^{t+1 }_{i,j} - S_y B_{i,j}C^{t+1 }_{i,j-1} = \\
|
||||
S_x A_{i,j} C^{t+1/2}_{i+1,j} + (1 - 2S_x\alpha^x_{i,j})C^{t+1/2}_{i,j} + S_x A_{i,j}C^{t+1/2}_{i-1,j}
|
||||
\end{split}
|
||||
\end{equation}
|
||||
|
||||
\noindent If the diffusion coefficients are constant,
|
||||
$A_{i,j}=B_{i,j}=\alpha$ and the scheme reverts to the homogeneous
|
||||
case. Problem with this discretization is that the terms in $A_{ij}$
|
||||
and $B_{ij}$ can be negative depending on the derivative of the
|
||||
diffusion coefficient, resulting in unphysical values for the
|
||||
concentrations.
|
||||
|
||||
** Direct discretization
|
||||
|
||||
As noted in literature (LeVeque and Numerical Recipes) a better way is
|
||||
to discretize directly the physical problem ([[eqn:hetdiff]]) at points
|
||||
halfway between grid points:
|
||||
|
||||
\begin{align*}
|
||||
\begin{cases}
|
||||
\displaystyle \alpha(x_{i+1/2}) \frac{\partial C }{\partial x}(x_{i+1/2}) & \displaystyle = \alpha_{i+1/2} \left( \frac{C_{i+1} -C_{i}}{\Delta x} \right) \\
|
||||
\displaystyle \alpha(x_{i-1/2}) \frac{\partial C }{\partial x}(x_{i-1/2}) & \displaystyle = \alpha_{i-1/2} \left( \frac{C_{i} -C_{i-1}}{\Delta x} \right)
|
||||
\end{cases}
|
||||
\end{align*}
|
||||
|
||||
\noindent A further differentiation gives us the centered
|
||||
approximation of $\frac{\partial}{\partial x} \left(\alpha(x)
|
||||
\frac{\partial C }{\partial x}\right)$:
|
||||
|
||||
\begin{align*}
|
||||
\frac{\partial}{\partial x} \left(\alpha(x)
|
||||
\frac{\partial C }{\partial x}\right)(x_i) & \simeq \frac{1}{\Delta x}\left[\alpha_{i+1/2} \left( \frac{C_{i+1} -C_{i}}{\Delta x} \right) - \alpha_{i-1/2} \left( \frac{C_{i} -C_{i-1}}{\Delta x} \right) \right]\\
|
||||
&\displaystyle =\frac{1}{\Delta x^2} \left[ \alpha_{i+1/2}C_{i+1} - (\alpha_{i+1/2}+\alpha_{i-1/2}) C_{i} + \alpha_{i-1/2}C_{i-1}\right]
|
||||
\end{align*}
|
||||
|
||||
\noindent The ADI scheme with this approach becomes:
|
||||
|
||||
#+NAME: eqn:genADI_hetdir
|
||||
\begin{equation}
|
||||
\left\{
|
||||
\begin{aligned}
|
||||
\frac{C^{t+1/2}_{i,j}-C^{t }_{i,j}}{\Delta t/2} = & \frac{1}{\Delta x^2} \left[ \alpha_{i+1/2,j} C^{t+1/2}_{i+1,j} - (\alpha_{i+1/2,j}+\alpha_{i-1/2,j}) C^{t+1/2}_{i,j} + \alpha_{i-1/2,j} C^{t+1/2}_{i-1,j}\right] + \\
|
||||
& \frac{1}{\Delta y^2} \left[ \alpha_{i,j+1/2}C^{t}_{i,j+1} - (\alpha_{i,j+1/2}+\alpha_{i,j-1/2}) C^t_{i,j} + \alpha_{i,j-1/2}C^{t}_{i,j-1}\right]\\
|
||||
\frac{C^{t+1 }_{i,j}-C^{t+1/2}_{i,j}}{\Delta t/2} = & \frac{1}{\Delta y^2} \left[ \alpha_{i+1/2,j}C^{t+1/2}_{i+1,j} - (\alpha_{i+1/2,j}+\alpha_{i-1/2,j}) C^t_{i,j} + \alpha_{i-1/2,j}C^{t+1/2}_{i-1,j}\right] + \\
|
||||
& \frac{1}{\Delta x^2} \left[ \alpha_{i,j+1/2}C^{t+1}_{i,j+1} - (\alpha_{i,j+1/2}+\alpha_{i,j-1/2}) C^{t+1}_{i,j} + \alpha_{i,j-1/2}C^{t+1}_{i,j-1}\right]
|
||||
\end{aligned}
|
||||
\right.
|
||||
\end{equation}
|
||||
|
||||
\noindent Doing the usual algebra and separating implicit from
|
||||
explicit terms, the two sweeps become:
|
||||
|
||||
#+NAME: eqn:sweepX_hetdir
|
||||
\begin{equation}
|
||||
\left\{
|
||||
\begin{aligned}
|
||||
-S_x \alpha^x_{i+1/2,j} C^{t+1/2}_{i+1,j} + (1 + S_x(\alpha^x_{i+1/2,j}+ \alpha^x_{i-1/2,j}))C^{t+1/2}_{i,j} - S_x \alpha^x_{i-1/2,j} C^{t+1/2}_{i-1,j} = \quad & \\
|
||||
S_y \alpha^y_{i,j+1/2} C^{t }_{i,j+1} + (1 - S_y(\alpha^y_{i,j+1/2}+ \alpha^y_{i,j-1/2}))C^{t }_{i,j} + S_y \alpha^y_{i,j-1/2} C^{t }_{i,j-1} & \\[1em]
|
||||
-S_y \alpha^y_{i,j+1/2} C^{t+1 }_{i,j+1} + (1 + S_y(\alpha^y_{i,j+1/2}+ \alpha^y_{i,j-1/2}))C^{t+1 }_{i,j} - S_y \alpha^y_{i,j-1/2} C^{t+1 }_{i,j-1} = \qquad & \\
|
||||
S_x \alpha^x_{i+1/2,j} C^{t+1/2}_{i+1,j} + (1 - S_x(\alpha^x_{i+1/2,j}+ \alpha^x_{i-1/2,j}))C^{t+1/2}_{i,j} + S_x \alpha^x_{i-1/2,j} C^{t+1/2}_{i-1,j}
|
||||
\end{aligned}
|
||||
\right.
|
||||
\end{equation}
|
||||
|
||||
@ -1,3 +1,4 @@
|
||||
## Time-stamp: "Last modified 2023-01-05 17:52:55 delucia"
|
||||
|
||||
## Brutal implementation of 2D ADI scheme
|
||||
## Square NxN grid with dx=dy=1
|
||||
@ -23,12 +24,12 @@ ADI <- function(n, dt, iter, alpha) {
|
||||
tmpY[i,] <- SweepByRow(i, resY, dt=dt, alpha=alpha)
|
||||
|
||||
res <- t(tmpY)
|
||||
out[[it]] <- res }
|
||||
out[[it]] <- res
|
||||
}
|
||||
|
||||
return(out)
|
||||
}
|
||||
|
||||
|
||||
## Workhorse function to fill A, B and solve for a given *row* of the
|
||||
## grid matrix
|
||||
SweepByRow <- function(i, field, dt, alpha) {
|
||||
@ -48,12 +49,12 @@ SweepByRow <- function(i, field, dt, alpha) {
|
||||
B <- numeric(ncol(field))
|
||||
|
||||
## We now distinguish the top and bottom rows
|
||||
if (i == 1){
|
||||
if (i == 1) {
|
||||
## top boundary, "i-1" doesn't exist or is at a ghost
|
||||
## node/cell boundary (TODO)
|
||||
for (ii in seq_along(B))
|
||||
B[ii] <- (-1 +2*Sy)*field[i,ii] - Sy*field[i+1,ii]
|
||||
} else if (i == nrow(field)){
|
||||
} else if (i == nrow(field)) {
|
||||
## bottom boundary, "i+1" doesn't exist or is at a ghost
|
||||
## node/cell boundary (TODO)
|
||||
for (ii in seq_along(B))
|
||||
@ -122,3 +123,240 @@ plot(adi2[[length(adi2)]], ref2, log="xy", xlab="ADI", ylab="ode.2D (reference)"
|
||||
las=1, xlim=c(1E-15, 1), ylim=c(1E-15, 1))
|
||||
abline(0,1)
|
||||
|
||||
|
||||
## Test heterogeneous scheme, chain rule
|
||||
ADIHet <- function(field, dt, iter, alpha) {
|
||||
|
||||
if (!all.equal(dim(field), dim(alpha)))
|
||||
stop("field and alpha are not matrix")
|
||||
|
||||
## now both field and alpha must be nx*ny matrices
|
||||
nx <- ncol(field)
|
||||
ny <- nrow(field)
|
||||
dx <- dy <- 1
|
||||
|
||||
## find out the center of the grid to apply conc=1
|
||||
cenx <- ceiling(nx/2)
|
||||
ceny <- ceiling(ny/2)
|
||||
field[cenx, ceny] <- 1
|
||||
|
||||
Aij <- Bij <- alpha
|
||||
|
||||
for (i in seq(2,ncol(field)-1)) {
|
||||
for (j in seq(2,nrow(field)-1)) {
|
||||
Aij[i,j] <- (alpha[i+1,j]-alpha[i-1,j])/4 + alpha[i,j]
|
||||
Bij[i,j] <- (alpha[i,j+1]-alpha[i,j-1])/4 + alpha[i,j]
|
||||
}
|
||||
}
|
||||
|
||||
if (any(Aij<0) || any(Bij<0))
|
||||
stop("Aij or Bij are negative!")
|
||||
|
||||
## prepare containers for computations and outputs
|
||||
tmpX <- tmpY <- res <- field
|
||||
out <- vector(mode="list", length=iter)
|
||||
|
||||
for (it in seq(1, iter)) {
|
||||
for (i in seq(1, ny))
|
||||
tmpX[i,] <- SweepByRowHet(i, res, dt=dt, alpha=alpha, Aij, Bij)
|
||||
|
||||
resY <- t(tmpX)
|
||||
for (i in seq(1, nx))
|
||||
tmpY[i,] <- SweepByRowHet(i, resY, dt=dt, alpha=alpha, Bij, Aij)
|
||||
|
||||
res <- t(tmpY)
|
||||
out[[it]] <- res
|
||||
}
|
||||
|
||||
return(out)
|
||||
}
|
||||
|
||||
|
||||
## Workhorse function to fill A, B and solve for a given *row* of the
|
||||
## grid matrix
|
||||
SweepByRowHet <- function(i, field, dt, alpha, Aij, Bij) {
|
||||
dx <- 1 ## fixed in our test
|
||||
Sx <- Sy <- dt/2/dx/dx
|
||||
|
||||
## diagonal of A at once
|
||||
A <- matrix(0, nrow(field), ncol(field))
|
||||
diag(A) <- 1+2*Sx*diag(alpha)
|
||||
|
||||
## adjacent diagonals "Sx"
|
||||
for (ii in seq(1, nrow(field)-1)) {
|
||||
A[ii+1, ii] <- -Sx*Aij[ii+1,ii]
|
||||
A[ii, ii+1] <- -Sx*Aij[ii,ii+1]
|
||||
}
|
||||
|
||||
B <- numeric(ncol(field))
|
||||
|
||||
## We now distinguish the top and bottom rows
|
||||
if (i == 1) {
|
||||
## top boundary, "i-1" doesn't exist or is at a ghost
|
||||
## node/cell boundary (TODO)
|
||||
for (ii in seq_along(B))
|
||||
B[ii] <- Sy*Bij[i+1,ii]*field[i+1,ii] + (1-2*Sy*Bij[i,ii])*field[i, ii]
|
||||
} else if (i == nrow(field)) {
|
||||
## bottom boundary, "i+1" doesn't exist or is at a ghost
|
||||
## node/cell boundary (TODO)
|
||||
for (ii in seq_along(B))
|
||||
B[ii] <- (1-2*Sy*Bij[i,ii])*field[i, ii] + Sy*Bij[i-1,ii]*field[i-1,ii]
|
||||
|
||||
} else {
|
||||
## inner grid row, full expression
|
||||
for (ii in seq_along(B))
|
||||
B[ii] <- Sy*Bij[i+1,ii]*field[i+1,ii] + (1-2*Sy*Bij[i,ii])*field[i, ii] + Sy*Bij[i-1,ii]*field[i-1,ii]
|
||||
}
|
||||
|
||||
x <- solve(A, B)
|
||||
x
|
||||
}
|
||||
|
||||
## adi2 <- ADI(n=51, dt=10, iter=200, alpha=1E-3)
|
||||
## ref2 <- DoRef(n=51, alpha=1E-3, dt=10, iter=200)
|
||||
|
||||
n <- 51
|
||||
field <- matrix(0, n, n)
|
||||
alphas <- matrix(1E-3*runif(n*n, 1,1.2), n, n)
|
||||
|
||||
## for (i in seq(1,nrow(alphas)))
|
||||
## alphas[i,] <- seq(1E-7,1E-3, length=n)
|
||||
|
||||
#diag(alphas) <- rep(1E-2, n)
|
||||
|
||||
adih1 <- ADIHet(field=field, dt=10, iter=100, alpha=alphas)
|
||||
adi2 <- ADI(n=n, dt=10, iter=100, alpha=1E-3)
|
||||
|
||||
|
||||
par(mfrow=c(1,3))
|
||||
image(adi2[[length(adi2)]])
|
||||
image(adih1[[length(adih1)]])
|
||||
points(0.5,0.5, col="red",pch=4)
|
||||
plot(adih1[[length(adih1)]], adi2[[length(adi2)]], pch=4, log="xy")
|
||||
abline(0,1)
|
||||
|
||||
|
||||
sapply(adih1, sum)
|
||||
sapply(adi2, sum)
|
||||
|
||||
adi2
|
||||
|
||||
|
||||
par(mfrow=c(1,2))
|
||||
image(alphas)
|
||||
image(adih1[[length(adih1)]])
|
||||
points(0.5,0.5, col="red",pch=4)
|
||||
|
||||
|
||||
|
||||
|
||||
## Test heterogeneous scheme, direct discretization
|
||||
ADIHetDir <- function(field, dt, iter, alpha) {
|
||||
|
||||
if (!all.equal(dim(field), dim(alpha)))
|
||||
stop("field and alpha are not matrix")
|
||||
|
||||
## now both field and alpha must be nx*ny matrices
|
||||
nx <- ncol(field)
|
||||
ny <- nrow(field)
|
||||
dx <- dy <- 1
|
||||
|
||||
## find out the center of the grid to apply conc=1
|
||||
cenx <- ceiling(nx/2)
|
||||
ceny <- ceiling(ny/2)
|
||||
field[cenx, ceny] <- 1
|
||||
|
||||
## prepare containers for computations and outputs
|
||||
tmpX <- tmpY <- res <- field
|
||||
out <- vector(mode="list", length=iter)
|
||||
|
||||
for (it in seq(1, iter)) {
|
||||
for (i in seq(2, ny-1)) {
|
||||
Aij <- cbind(harm(alpha[i,], alpha[i-1,]), harm(alpha[i,], alpha[i+1,]))
|
||||
Bij <- cbind(harm(alpha[,i], alpha[,i-1]), harm(alpha[,i], alpha[,i+1]))
|
||||
tmpX[i,] <- SweepByRowHetDir(i, res, dt=dt, Aij, Bij)
|
||||
}
|
||||
resY <- t(tmpX)
|
||||
for (i in seq(2, nx-1))
|
||||
tmpY[i,] <- SweepByRowHetDir(i, resY, dt=dt, Bij, Aij)
|
||||
res <- t(tmpY)
|
||||
out[[it]] <- res
|
||||
}
|
||||
|
||||
return(out)
|
||||
}
|
||||
|
||||
harm <- function(x,y) 1/(1/x+1/y)
|
||||
|
||||
harm(1,4)
|
||||
|
||||
|
||||
## Direct discretization, Workhorse function to fill A, B and solve
|
||||
## for a given *row* of the grid matrix
|
||||
SweepByRowHetDir <- function(i, field, dt, Aij, Bij) {
|
||||
dx <- 1 ## fixed in our test
|
||||
Sx <- Sy <- dt/2/dx/dx
|
||||
|
||||
## diagonal of A at once
|
||||
A <- matrix(0, nrow(field), ncol(field))
|
||||
diag(A) <- 1 + Sx*(Aij[,1]+Aij[,2])
|
||||
|
||||
## adjacent diagonals "Sx"
|
||||
for (ii in seq(1, nrow(field)-1)) {
|
||||
A[ii+1, ii] <- -Sx*Aij[ii,2] # i-1/2
|
||||
A[ii, ii+1] <- -Sx*Aij[ii,1] # i+1/2
|
||||
}
|
||||
|
||||
B <- numeric(ncol(field))
|
||||
|
||||
for (ii in seq_along(B))
|
||||
B[ii] <- Sy*Bij[ii,2]*field[i+1,ii] + (1 - Sy*(Bij[ii,1]+Bij[ii,2]))*field[i, ii] + Sy*Bij[ii,1]*field[i-1,ii]
|
||||
|
||||
lastA <<- A
|
||||
lastB <<- B
|
||||
x <- solve(A, B)
|
||||
x
|
||||
}
|
||||
|
||||
## adi2 <- ADI(n=51, dt=10, iter=200, alpha=1E-3)
|
||||
## ref2 <- DoRef(n=51, alpha=1E-3, dt=10, iter=200)
|
||||
|
||||
n <- 5
|
||||
field <- matrix(0, n, n)
|
||||
alphas <- matrix(1E-5*runif(n*n, 1,2), n, n)
|
||||
|
||||
alphas1 <- matrix(3E-5, n, 25)
|
||||
alphas2 <- matrix(1E-5, n, 26)
|
||||
|
||||
alphas <- cbind(alphas1, alphas2)
|
||||
|
||||
## for (i in seq(1,nrow(alphas)))
|
||||
## alphas[i,] <- seq(1E-7,1E-3, length=n)
|
||||
|
||||
#diag(alphas) <- rep(1E-2, n)
|
||||
|
||||
adih <- ADIHetDir(field=field, dt=10, iter=200, alpha=alphas)
|
||||
adi2 <- ADI(n=n, dt=10, iter=200, alpha=1E-5)
|
||||
|
||||
|
||||
par(mfrow=c(1,3))
|
||||
image(adi2[[length(adi2)]])
|
||||
image(adih[[length(adih)]])
|
||||
points(0.5,0.5, col="red",pch=4)
|
||||
plot(adih[[length(adih)]], adi2[[length(adi2)]], pch=4, log="xy")
|
||||
abline(0,1)
|
||||
|
||||
|
||||
sapply(adih, sum)
|
||||
sapply(adi2, sum)
|
||||
|
||||
adi2
|
||||
|
||||
|
||||
par(mfrow=c(1,2))
|
||||
image(alphas)
|
||||
points(0.5,0.5, col="red",pch=4)
|
||||
image(adih[[length(adih)]])
|
||||
points(0.5,0.5, col="red",pch=4)
|
||||
|
||||
options(width=110)
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user