Gruppeoppgaver for uke 8

Tilbake til hovedsiden
Tilbake til oversikten over l?sningsforslag

F?lgende R-kode l?ser differenslikningen i oppgave 3.12.


#Steffen Gr?nneberg 20.02.06 - steffeng@math.uio.no

N <- 20
beta <- 0.01
b <- 0.0025
gamma <- 0.0025


D <- 0*(0:N) %*% t(0:N)
PI <- beta* ((N*seq(1,N-2) - (seq(1,N-2))^2)/N)

# First we generate the subdiagonal elements
X <- -seq(1,N-2)*(b + gamma)
X <- c(X, -N*(b+gamma))
X <- cbind(rbind(0,diag(X)),0)

# Then the superdiagonal
Z <- c(0,-PI)
Z <- cbind(0,rbind(diag(Z),0))

# And finally the diagonal ones.
Y <- PI + ((seq(1,N-2))*(b + gamma))
Y <- c(1,Y, N*(b+gamma))
Y <- diag(Y)

D <- X + Y + Z

cvec <- c(0, 1-PI,1)

tau <- solve(D) %*% cvec

plot(tau, type='l', xlab='k' ,main='Forventet tid til absorbsjon n?r man starter p? k')

Merk at R lager nye kolonner eller rader hvis det trengs for cbind og rbind operasjonene!

Dessuten, solve(D) gir den inverse av D.