The question is solved with R and this report is made by R Markdown. One of the reasons is that they are both in RStudio, which is new to me and I want to try. In R Markdown, the Chinese characters are so ugly and glare so I decided to use English to convey what I want to express. This is my first report written in English and shaped by R with R Markdown!

Explicit Difference Method

Flowing is the code chunk based on the explicit difference method, most codes of which are in a function named calex and you can call calex function to solve the problems with alternative different parameters. And you can change the parameters as you want.

I use backward difference method because it is easily for codeing, comparing to central difference method which is mainly discussed in class. What you should do is just adding one column named x = -1 in the initial data frame.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# Backward Difference Method
calex <- function(dx,dt,BOD,x = 8,u = 5,E = 2, k1 = 0.0151,t = 1,plot = FALSE) {
  col <- x/dx+3; row <- t/dt+1; name <- c('t',seq(-1,x,dx))
  alpha <- E*dt/dx^2; beta <- u*dt/dx-2*dt*E/dx^2-k1*dt; gamma <- 1-u*dt/dx+E*dt/dx^2
  ex <- rep(0,row*col);attr(ex,"dim") <- c(row,col); ex <- as.tibble(ex)
  ex[,1] <- seq(0.0,t,dt); ex[,2] <- BOD; ex[,3] <- BOD
  if (dt*u/dx <= 1 && E*dt/dx^2<= 0.5) {
    for (j in 1:(t/dt)) for (i in 1:(x/dx)) 
      ex[j+1,i+3] <- alpha*ex[j,i+1]+beta*ex[j,i+2]+gamma*ex[j,i+3]
    if (plot) {
      names(ex) <- name
      pex <- ex %>%
        dplyr::filter(t != 0) %>%
        melt(id = "t",variable.name = "x",value.name = "BOD") %>%
        dplyr::filter(x != -1) %>%
        ggvis(~x,~t) %>%
        layer_points(size:=~factor(BOD*25))
      return(list(ex,pex))
    } else {
      names(ex) <- name
      return(ex)
    }
  } else return(cat("You should choose proper values for both dx and dt !"))
}

And you can call this function with the specific paremeters:
This is what the question requires:

1
2
3
4
# This is what the question requires.
knitr::kable(calex(1,0.2,20),"html",align = "c",table.attr = "class=\"table table-bordered\"") %>%
  kable_styling(bootstrap_options = c("striped","hover","responsive")) %>%
  add_header_above(c("","x" = length(calex(1,0.2,20))-1))
1
2
## Warning: `as.tibble()` is deprecated, use `as_tibble()` (but mind the new semantics).
## This warning is displayed once per session.
x
t -1 0 1 2 3 4 5 6 7 8
0.0 20 20 0.00000 0.00000 0.00000 0.000000 0.000000 0.000000 0.000000 0.000000
0.2 20 20 11.93960 8.00000 0.00000 0.000000 0.000000 0.000000 0.000000 0.000000
0.4 20 20 16.71544 13.55186 6.35168 3.200000 0.000000 0.000000 0.000000 0.000000
0.6 20 20 18.62578 16.71335 11.89629 7.951899 3.171008 1.280000 0.000000 0.000000
0.8 20 20 19.38991 18.35425 15.50102 12.209433 7.593286 4.317385 1.520538 0.512000
1.0 20 20 19.69556 19.16112 17.57179 15.278863 11.642738 8.106452 4.495968 2.231269

In this function as we can see, 1 is the quantity of dx, 0.2 is the quantity of dt and 20 is the quantity of BOD. And in this table, the numeric row names means the value of x(Km). Pay attention to the -1 column, this column is added for backward differece method.

You can change them as you want:

1
2
3
knitr::kable(calex(1,0.2,10),table.attr = "class=\"table table-bordered\"") %>%
  kable_styling(bootstrap_options = c("striped","hover","responsive")) %>%
  add_header_above(c("","x" = length(calex(1,0.2,20))-1))
x
t -1 0 1 2 3 4 5 6 7 8
0.0 10 10 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.0000000 0.000000
0.2 10 10 5.969800 4.000000 0.000000 0.000000 0.000000 0.000000 0.0000000 0.000000
0.4 10 10 8.357720 6.775931 3.175840 1.600000 0.000000 0.000000 0.0000000 0.000000
0.6 10 10 9.312888 8.356676 5.948147 3.975949 1.585504 0.640000 0.0000000 0.000000
0.8 10 10 9.694955 9.177123 7.750512 6.104716 3.796643 2.158692 0.7602688 0.256000
1.0 10 10 9.847782 9.580561 8.785897 7.639432 5.821369 4.053226 2.2479839 1.115635

Sometimes you want to change the default value of, for example, x , k1 or t. You can do it because they are all parameters, which have default values perserved in the function beforehand. Just like this:

1
2
knitr::kable(calex(0.5,0.1,10,t = 2,x = 10),align = "c",
             table.attr = "class=\"table table-bordered\"")
1
## You should choose proper values for both dx and dt !

Error! But don’t warry! Your values of dx and dt just fail to meet the condition. Change them like this:

And this is the same table as the corresponding example in our book:

1
2
3
4
knitr::kable(calex(1,0.2,10,t = 2,x = 10),"html",align = "c",
             table.attr = "class=\"table table-bordered\"") %>%
  kable_styling(bootstrap_options = c("striped","hover","responsive")) %>%
  add_header_above(c("","x" = length(calex(1,0.2,10,t = 2,x = 10))-1))
x
t -1 0 1 2 3 4 5 6 7 8 9 10
0.0 10 10 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.0000000 0.000000 0.0000000 0.0000000
0.2 10 10 5.969800 4.000000 0.000000 0.000000 0.000000 0.000000 0.0000000 0.000000 0.0000000 0.0000000
0.4 10 10 8.357720 6.775931 3.175840 1.600000 0.000000 0.000000 0.0000000 0.000000 0.0000000 0.0000000
0.6 10 10 9.312888 8.356676 5.948147 3.975949 1.585504 0.640000 0.0000000 0.000000 0.0000000 0.0000000
0.8 10 10 9.694955 9.177123 7.750512 6.104716 3.796643 2.158692 0.7602688 0.256000 0.0000000 0.0000000
1.0 10 10 9.847782 9.580561 8.785897 7.639432 5.821369 4.053226 2.2479839 1.115635 0.3545344 0.1024000
1.2 10 10 9.908913 9.772041 9.340651 8.618643 7.347721 5.823756 4.0261456 2.510352 1.2607650 0.5570501
1.4 10 10 9.933365 9.860674 9.624722 9.196195 8.373049 7.224314 5.6967104 4.126714 2.6092534 1.4753064
1.6 10 10 9.943146 9.900944 9.765590 9.518625 9.010575 8.217527 7.0509492 5.662549 4.1352656 2.7547787
1.8 10 10 9.947058 9.918978 9.833782 9.691454 9.385445 8.869364 8.0432980 6.940926 5.5898948 4.1814957
2.0 10 10 9.948623 9.926963 9.866177 9.781231 9.596713 9.273072 8.7185845 7.908485 6.8205008 5.5500663

Besides, the visualization of the result can be shown here:

1
calex(1,0.2,10,t = 2,x = 10,plot = TRUE)[[2]]

As you can see, the size of each point means the quantity of the BOD value. And I ignore the condition where x = -1 and t = 0.

Implicit Difference Method

Next comes the implicit differece method based on calim function. You can call this function with different parameters to solve the specific problems. As you can see, the method is more complex and difficult than the last method.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
calim <- function(dx,dt,BOD,x = 8,u = 5,E = 2, k1 = 0.0151,t = 1,plot = FALSE) {
  col <- x/dx+2; row <- t/dt+1; name <- c('t',seq(0,x,dx))
  alpha <- -E/dx^2; beta <- 1/dt+2*E/dx^2+k1/2; gamma <- alpha
  d2 <- 1/dt-u/dx; d1 <- u/dx-k1/2
  im <- rep(0,row*col);attr(im,"dim") <- c(row,col); im <- as.tibble(im)
  im[,1] <- seq(0.0,t,dt); im[,2] <- BOD
  delta <- rep(0,(row-1)*(col-1-1)); attr(delta,'dim') <- c(row-1,col-1-1); delta <- as.tibble(delta)
  g <- delta[1:(row-1),1:(col-1-1)]; omega <- delta[1:(row-1),1:(col-1-1-1)]
  if(dt*u/dx > 1 | dt*u/dx <= 0.5) {
    return(cat("You should choose proper values for both dx and dt !"))
  } else {
    for(j in 1:(t/dt)) {
      delta[j,1] <- im[j,2]*d1+im[j,3]*d2
      delta[j,1] <- delta[j,1]-alpha*im[j,2]
      delta[j,2:(x/dx)] <- im[j,3:(x/dx+1)]*d1+im[j,4:(x/dx+2)]*d2
      for(i in 1:(x/dx-2)) {
        omega[j,1] <- gamma/beta; g[j,1] <- delta[j,1]/beta
        g[j,i+1] <- (delta[j,i+1]-alpha*g[j,i])/(beta-alpha*omega[j,i])
        omega[j,i+1] <- gamma/(beta-alpha*omega[j,i])
      }
      alphap <- alpha - gamma; betap <- beta + 2*gamma
      g[j,x/dx] <- (delta[j,x/dx]-alphap*g[j,x/dx-1])/(betap-alphap*omega[j,x/dx-1])
      im[j+1,x/dx+2] <- g[j,x/dx]
      for(k in (x/dx+1):3) im[j+1,k] <- -im[j+1,k+1]*omega[j,k-2]+g[j,k-2]
    }
    if (plot) {
      names(im) <- name
      pim <- im %>%
        dplyr::filter(t != 0) %>%
        melt(id = "t",variable.name = "x",value.name = "BOD") %>%
        dplyr::filter(x != 0) %>%
        ggvis(~x,~t) %>%
        layer_points(size:=~factor(BOD*25))
      return(list(im,pim))
    } else {
      names(im) <- name
      return(im)
    }
  }
}

And you can call this function with the specific paremeters:
This is what the question requires:

1
2
3
4
5
knitr::kable(calim(dx = 0.5,dt = 0.1, BOD = 20),"html",align = "c",
             table.attr = "class=\"table table-bordered\"") %>%
  kable_styling(bootstrap_options = c("striped","hover","responsive")) %>%
  add_header_above(c("","x" = length(calim(dx = 0.5,dt = 0.1, BOD = 20))-1)) %>%
  scroll_box(width = "100%")
x
t 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8
0.0 20 0.00000 0.000000 0.000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
0.1 20 15.47371 5.323023 1.831144 0.6299215 0.2166958 0.0745443 0.0256436 0.0088215 0.0030346 0.0010439 0.0003591 0.0001235 0.0000424 0.0000144 0.0000044 0.0000000
0.2 20 18.06790 13.756609 7.326530 3.4127739 1.4810040 0.6150795 0.2479195 0.0977829 0.0379369 0.0145294 0.0055068 0.0020689 0.0007707 0.0002822 0.0000937 0.0000044
0.3 20 18.99612 16.774201 12.968034 8.2013711 4.5429130 2.3046370 1.0994749 0.5014252 0.2209647 0.0947823 0.0397817 0.0163979 0.0066488 0.0026326 0.0009470 0.0000936
0.4 20 19.42205 18.158863 15.884167 12.5277515 8.6450437 5.3328253 3.0173174 1.5976826 0.8033530 0.3876644 0.1809250 0.0821242 0.0363666 0.0156195 0.0061069 0.0009456
0.5 20 19.64369 18.879391 17.472921 15.2426522 12.2398650 8.9005941 5.8973198 3.6102622 2.0706404 1.1256777 0.5854414 0.2933448 0.1422205 0.0664284 0.0283105 0.0060977
0.6 20 19.76811 19.283897 18.386685 16.9087607 14.7580832 12.0300076 9.0625149 6.3143759 4.0990819 2.5020880 1.4487193 0.8015823 0.4259303 0.2166890 0.1008721 0.0282678
0.7 20 19.84160 19.522808 18.934465 17.9453990 16.4390139 14.3769174 11.8658597 9.1721689 6.6327546 4.5035325 2.8879926 1.7599202 1.0238769 0.5674248 0.2887781 0.1007199
0.8 20 19.88661 19.669128 19.273340 18.6023308 17.5515793 16.0420534 14.0669809 11.7313318 9.2497928 6.8826694 4.8406946 3.2289910 2.0493083 1.2349557 0.6865832 0.2883424
0.9 20 19.91494 19.761212 19.488158 19.0258468 18.2903135 17.1995760 15.7016309 13.8081156 11.6173257 9.3060635 7.0826517 5.1223975 3.5236746 2.2996790 1.3927511 0.6855473
1.0 20 19.93314 19.820404 19.626999 19.3030010 18.7841651 17.9989104 16.8836498 15.4056475 13.5870236 11.5178887 9.3462915 7.2425826 5.3523068 3.7593031 2.4677084 1.3906496

Similarly, you can change the parameters:
And this is the same table as the corresponding example in our book:

1
2
3
4
5
knitr::kable(calim(dx = 0.5,dt = 0.1, BOD = 10, x = 8),"html",align = "c",
             table.attr = "class=\"table table-bordered\"")%>%
  kable_styling(bootstrap_options = c("striped","hover","responsive")) %>%
  add_header_above(c("","x" = length(calim(dx = 0.5,dt = 0.1, BOD = 10, x = 8))-1)) %>%
  scroll_box(width = "100%")
x
t 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8
0.0 10 0.000000 0.000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
0.1 10 7.736853 2.661512 0.9155718 0.3149608 0.1083479 0.0372722 0.0128218 0.0044108 0.0015173 0.0005220 0.0001796 0.0000618 0.0000212 0.0000072 0.0000022 0.0000000
0.2 10 9.033951 6.878304 3.6632653 1.7063870 0.7405020 0.3075397 0.1239598 0.0488915 0.0189685 0.0072647 0.0027534 0.0010345 0.0003853 0.0001411 0.0000469 0.0000022
0.3 10 9.498061 8.387100 6.4840172 4.1006856 2.2714565 1.1523185 0.5497374 0.2507126 0.1104823 0.0473911 0.0198909 0.0081990 0.0033244 0.0013163 0.0004735 0.0000468
0.4 10 9.711024 9.079431 7.9420837 6.2638758 4.3225218 2.6664126 1.5086587 0.7988413 0.4016765 0.1938322 0.0904625 0.0410621 0.0181833 0.0078098 0.0030535 0.0004728
0.5 10 9.821843 9.439695 8.7364606 7.6213261 6.1199325 4.4502970 2.9486599 1.8051311 1.0353202 0.5628388 0.2927207 0.1466724 0.0711103 0.0332142 0.0141552 0.0030489
0.6 10 9.884056 9.641949 9.1933426 8.4543803 7.3790416 6.0150038 4.5312574 3.1571880 2.0495409 1.2510440 0.7243597 0.4007911 0.2129652 0.1083445 0.0504361 0.0141339
0.7 10 9.920801 9.761404 9.4672324 8.9726995 8.2195070 7.1884587 5.9329298 4.5860844 3.3163773 2.2517663 1.4439963 0.8799601 0.5119385 0.2837124 0.1443891 0.0503600
0.8 10 9.943305 9.834564 9.6366700 9.3011654 8.7757896 8.0210267 7.0334905 5.8656659 4.6248964 3.4413347 2.4203473 1.6144955 1.0246541 0.6174778 0.3432916 0.1441712
0.9 10 9.957468 9.880606 9.7440791 9.5129234 9.1451567 8.5997880 7.8508154 6.9040578 5.8086629 4.6530317 3.5413258 2.5611987 1.7618373 1.1498395 0.6963755 0.3427736
1.0 10 9.966572 9.910202 9.8134997 9.6515005 9.3920826 8.9994552 8.4418249 7.7028238 6.7935118 5.7589444 4.6731458 3.6212913 2.6761534 1.8796515 1.2338542 0.6953248

Sometimes your values are not suitable:

1
2
knitr::kable(calim(dx = 1,dt = 0.1, BOD = 20),"html",align = "c",
             table.attr = "class=\"table table-bordered\"")
1
## You should choose proper values for both dx and dt !

And you can modify them until they are proper.

The visualization of the result can be available in the same way:

1
calim(dx = 0.5,dt = 0.1, BOD = 20,plot = TRUE)[[2]]

The condition where x = 0 and t = 0 are ignored.

Difference

So what’s the differences between .Rmd file and .Rmarkdown file?

  • HTML widgets are not supported.

  • You cannot use Markdown features only supported by Pandoc, such as citations. Math expressions only work if you have installed the xaringan package (Xie 2018e) and applied the JavaScript solution mentioned in Section B.3.(See more)

  • The tabulate in .Rmarkdown file is similar to .md file, and toc can be generated in .Rmarkdown file.

  • The color of code block in .Rmarkdown file is similar to .md file.