Machanie

Funkcja do interpolacji (smoothing spline):

approxfun2 = function(x,y,...) {
  m = smooth.spline(x,y,...)
  function(xnew) {
    ynew = predict(m, xnew)$y
    ynew[xnew>max(x)] = NA
    ynew[xnew<min(x)] = NA
    ynew
  }
}

Ładujemy dane z Fluenta i robimy z nich funkcje czasu:

tab=read.table("~/cl",skip=1,header=TRUE)
tab = tab[-(1:37),]
cl_mr = approxfun2(tab$Flow.Time,tab$cl)

Ładujemy dane z TCLB i robimy z nich funkcje czasu:

tab=read.csv("~/long_base__Log_P00_00000000.csv")
f = 1/1500
R = tab$PDX*2
tab$umax = 2*pi*f*R
tab$cl = -2*tab$ForceY/(tab$PDX * tab$umax^2)
cl_ll = approxfun2(tab$Time_si,tab$cl)

Funkcja do analizy:

decompose = function(fun, k = 1500/2, dt = 5e-2/2) {
  A = sapply(0:30,function(i) {
    t = seq(0,dt,len=k+1)[-1]+dt*i
    fun(t)
  }) # Dzieli wyniki na kawalki
  sel = !apply(is.na(A),2,any) # eliminuje kawalki z poza zakresu
  A = A[,sel]
  L = colMeans(A) # Liczy Lrednie w kazdym kawalku
  O = rowMeans(A) # Liczy Lrednie po kawalkach
  O.sd = apply(A,1,sd) # Liczy odchylenie standardowe po kawalkach
  A = A - O # Liczy reszta
  pca = prcomp(t(A)) # liczy PCA/POD tej reszty.
  var = pca$sdev^2
  list(
    n = ncol(A),
    cuts = A,
    means = L,
    trend = O,
    trend.sd = O.sd,
    pc.var = var,
    pc.shape = pca$rotation,
    pc.coef = t(A) %*% pca$rotation,
    pc90 = sum(cumsum(var)/sum(var) < 0.90)+1,
    pc99 = sum(cumsum(var)/sum(var) < 0.99)+1
  )
}    

Puszczenie analizy na obu zestawach danych:

dec_mr = decompose(cl_mr)
dec_ll = decompose(cl_ll)

Plot średniego \(C_L\) w pół-okresach:

plot(dec_mr$means, type="b", ylim=c(-1,3), xlab="Half period",ylab="CL")
lines(dec_ll$means, type="b", col=2)

Plot \(C_L\) w pół-okresie (z \(99\%\) przedziałem ufności):

p = seq(0,0.5,len=1500/2+1)[-1]
s = qnorm(0.995)
plot(p,dec_mr$trend, type="l", ylim=c(-1,3))
lines(p,dec_mr$trend + dec_mr$trend.sd*s, lty=2)
lines(p,dec_mr$trend - dec_mr$trend.sd*s, lty=2)
lines(p,dec_ll$trend, col=2)
lines(p,dec_ll$trend + dec_ll$trend.sd*s, lty=2, col=2)
lines(p,dec_ll$trend - dec_ll$trend.sd*s, lty=2, col=2)

Zróbmy analize PCA/POD reszt (czyli przebiegów minus średni przebieg):

ret = prcomp(t(cbind(dec_mr$cuts/dec_mr$n,dec_ll$cuts/dec_ll$n)))
col = gray.colors(length(ret$sdev))
col[1:3]= 1:3+1
barplot(ret$sdev,ylab="Standard Deviation",names.arg = colnames(ret$rotation),col=col)

I zobaczmy jaką część zmienności realizują kolejne mody (kształty z PCA/POD):

W = cbind(apply(t(dec_mr$cuts) %*% ret$rotation,2,var),
          apply(t(dec_ll$cuts) %*% ret$rotation,2,var))
barplot(W,col=col,names.arg = c("MR","LLW"), ylab="Total variance")

Widać że pierwsze 3 mody dominują. wyglądają one tak:

matplot(p,ret$rotation[,1:3],col=1:3+1,type="l",lty=1)