#-------------------------------# # week4: Abschluss Rechnen # # mit Vektoren und Matrizen # #-------------------------------# # Teil a+b hatten wir schon in dem # week3a+b.txt gemacht, nehmen wir # das als Vorlage: # a) eps = 2 x0 = 1 T = 50 dt = 0.01 tt = seq(from=0,to=T,by=dt) xt = x0*cos(eps*tt) plot(tt,xt) # b) T = 50 dt = 0.01 tt = seq(from=0,to=T,by=dt) M = length(tt) xt = rep(0,M) vt = rep(0,M) xt[1] = x0 vt[1] = 0 t0 = Sys.time() # wir messen die Rechenzeit for(k in 2:M) { xt[k] = xt[k-1] + vt[k-1]*dt vt[k] = vt[k-1] - eps^2*xt[k-1]*dt } Sys.time() - t0 xt_exact = x0*cos(eps*tt) plot(tt,xt) points(tt,xt_exact,col="red") # c) y = c(x0,0) xt = rep(0,M) xt[1] = x0 r1 = c(1,dt) r2 = c(-eps^2*dt,1) R = rbind(r1,r2) R for(k in 2:M) { y = R %*% y xt[k] = y[1] } xt_exact = x0*cos(eps*tt) plot(tt,xt) points(tt,xt_exact,col="red") # d) T = 50 dt = 0.01 tt = seq(from=0,to=T,by=dt) M = length(tt) A = rbind( c(0,1) , c(-eps^2,0) ) Id = diag(rep(1,2)) Id A2 = A %*% A A3 = A2 %*% A A4 = A2 %*% A2 expdtA2 = Id + dt*A + dt^2/2 * A2 expdtA3 = expdtA2 + dt^3/6 * A3 expdtA4 = expdtA3 + dt^4/24 * A4 y = c(x0,0) xt = rep(0,M) xt[1] = x0 for(k in 2:M) { y = expdtA3 %*% y # das ist also numerisch deutlich stabiler xt[k] = y[1] } xt_exact = x0*cos(eps*tt) plot(tt,xt) points(tt,xt_exact,col="red") # e) y = c(x0,0) xt = rep(0,M) xt[1] = x0 u1 = c( cos(eps*dt) , sin(eps*dt)/eps ) u2 = c( -eps*sin(eps*dt) , cos(eps*dt) ) U = rbind( u1 , u2 ) for(k in 2:M) { y = U %*% y xt[k] = y[1] } xt_exact = x0*cos(eps*tt) plot(tt,xt) points(tt,xt_exact,col="red")