利用R语言求解PDE中的热方程
heatfd <- function(xl, xr, yb, yt, M, N){
   f <- function(x) sin(2 * pi * x)^2
   l <- function(t) 0 * t
   r <- function(t) 0 * t
   D <- 1
   h <- (xr - xl)/M
   k <- (yt - yb)/N
   m <- M - 1
   n <- N
   sigma <- D * k/(h * h)
   a <- Diag(c(1 - 2 * sigma * ones(1,m)),0) + Diag(c(sigma * ones(m-1, 1)),1)
   a <- a + Diag(c(sigma * ones(m-1,1)), -1)
   lside <- l(yb + (0:n)*k)
   rside <- r(yb + (0:n)*k)
   w <- matrix(1,nrow=m,ncol=n+1)
   w[,1] <- f(xl + (1:m)*h)
   for (j in 1:n){
      w[,j+1] <- a %*% w[,j] + sigma*rbind(lside[j],zeros(m-2,1),rside[j])
   }
   w <- rbind(lside,w,rside)
   x <- (0:(m+1))*h
   t <- (0:n)*k
   xt <- mesh(x,t)    #library(plot3D)
   surf3D(xt$x,xt$y,w,colkey = F,border = "black",bty = "b2")
}

heatfd(1,0,1,0,10,250)