# pomůcka pro kreslení obrázku do 14. kapitoly # NLR.priklad <- function(soubor="../obr15_4.ps",theta0=-0.3,theta1=-.9,theta2=0.0, x1=2,x2=6,delka=0.15,h=1,npoints2=1, derivL=TRUE){ # funkční zadání, ale theta a beta podobná # ukázkový obrázek NLR # přepracováno 051215, 071219, 090106, 100103 # beta0 <- exp(theta0); beta1 <- exp(theta1); beta2<- exp(theta2) kriv <- krivost(theta0,x1=x1,x2=x2) krivT <- krivostT(beta0,x1=x1,x2=x2) npoints <- npoints2*2+1 c <- c(delka/sqrt(crossprod(c(fce1(theta0,x1=x1,x2=x2))))) #c <- c(h*delka/sqrt(crossprod(c(h*fce1(theta0,x1=x1,x2=x2))))) tt <- theta0+(-npoints2:npoints2)*c # cT <- c(delka/sqrt(crossprod(c(fceT1(beta0,x1=x1,x2=x2))))) bb <- beta0+(-npoints2:npoints2)*cT tt.lines<-seq(theta1,theta2,length=101) bb.lines<-seq(beta1,beta2,length=101) opar = par(mfrow=c(1,2)) # dve krivky par(mfg=c(1,1)) plot(fce0(tt.lines,x1=x1,x2=x2),asp=1,type="l", xlab=expression(y[1]),ylab=expression(y[2]),xlim=c(0.3,1.2), main=expression(theta)) par(mfg=c(1,2)) plot(fceT0(bb.lines,x1=x1,x2=x2),asp=1,type="l", xlab=expression(y[1]),ylab=expression(y[2]),xlim=c(0.3,1.2), main=expression(beta)) # tecne nadroviny par(mfg=c(1,1)) lines(matrix(fce0(theta0,x1=x1,x2=x2),npoints,2, byrow=T)+diag(tt-theta0)%*%matrix(h*fce1(theta0,x1=x1,x2=x2), npoints,2,byrow=T)) points(matrix(fce0(theta0,x1=x1,x2=x2),npoints,2, byrow=T)+diag(tt-theta0)%*%matrix(h*fce1(theta0,x1=x1,x2=x2), npoints,2,byrow=T),pch="+") par(mfg=c(1,2)) lines(matrix(fceT0(beta0,x1=x1,x2=x2),npoints,2,byrow=T) +diag(bb-beta0) %*%matrix(h*fceT1(beta0,x1=x1,x2=x2),npoints,2,byrow=T)) points(matrix(fceT0(beta0,x1=x1,x2=x2),npoints,2,byrow=T) +diag(bb-beta0) %*%matrix(h*fceT1(beta0,x1=x1,x2=x2),npoints,2,byrow=T),pch="+") par("lwd" = 2) # ještě naznačení derivací if (derivL){ par(mfg=c(1,1)) # bod dotyku eta0 <- c(fce0(theta0,x1=x1,x2=x2)) text(eta0[1],eta0[2],label="A",pos=3) # lineární aproximace eta00 <- c(fce0(theta0,x1=x1,x2=x2)+c*h*fce1(theta0,x1=x1,x2=x2)) text(eta00[1],eta00[2],lab="B",pos=1) # kvadratická aproximace eta000 <- eta00+c^2*h^2*fce2(theta0,x1=x1,x2=x2)/2 # bod na křivce text(fce0(theta0+h*c,x1=x1,x2=x2),lab="D",pos=2) points(fce0(theta0+h*c,x1=x1,x2=x2),pch="+") # kvadratická aproximace text(eta000,lab="C",pos=4) points(rbind(eta00,eta000),pch="+") lines(rbind(eta00,eta000),col="blue") # složka normálová text(eta00+c^2*h^2*fce2N(theta0,x1=x1,x2=x2)/2,lab="Cn",pos=2) points(rbind(eta00,eta00+c^2*h^2*fce2N(theta0,x1=x1,x2=x2)/2),pch="+") lines(rbind(eta00,eta00+c^2*h^2*fce2N(theta0,x1=x1,x2=x2)/2),col="red") # složka tečná text(eta00+c^2*h^2*fce2T(theta0,x1=x1,x2=x2)/2,lab="Ct",pos=4) points(rbind(eta00,eta00+c^2*h^2*fce2T(theta0,x1=x1,x2=x2)/2),pch="+") lines(rbind(eta00,eta00+c^2*h^2*fce2T(theta0,x1=x1,x2=x2)/2),col="cyan") cat("vzdálenost AB, param. theta", sqrt(crossprod(eta0-eta00)),"\n") par(mfg=c(1,2)) etaT0 <- c(fceT0(beta0,x1=x1,x2=x2)) text(etaT0[1],etaT0[2],label="A",pos=3) etaT00 <- c(fceT0(beta0,x1=x1,x2=x2)+cT*h*fceT1(beta0,x1=x1,x2=x2)) text(etaT00[1],etaT00[2],lab="B",pos=1) etaT000 <- etaT00+cT^2*h^2*fceT2(beta0,x1=x1,x2=x2)/2 text(etaT00+cT^2*h^2*fceT2(beta0,x1=x1,x2=x2)/2,lab="C",pos=4) text(etaT00+cT^2*h^2*fceT2N(beta0,x1=x1,x2=x2)/2,lab="Cn",pos=2) text(etaT00+cT^2*h^2*fceT2T(beta0,x1=x1,x2=x2)/2,lab="Ct",pos=4) text(fceT0(beta0+h*cT,x1=x1,x2=x2),lab="D",pos=2) points(fceT0(beta0+h*cT,x1=x1,x2=x2),pch="+") lines(rbind(etaT00,etaT00+cT^2*h^2*fceT2(beta0,x1=x1,x2=x2)/2),col="blue") lines(rbind(etaT00,etaT00+cT^2*h^2*fceT2N(beta0,x1=x1,x2=x2)/2),col="red") lines(rbind(etaT00,etaT00+cT^2*h^2*fceT2T(beta0,x1=x1,x2=x2)/2),col="cyan") points(rbind(etaT00,etaT00+cT^2*h^2*fceT2(beta0,x1=x1,x2=x2)/2),pch="+") points(rbind(etaT00,etaT00+cT^2*h^2*fceT2N(beta0,x1=x1,x2=x2)/2),pch="+") points(rbind(etaT00,etaT00+cT^2*h^2*fceT2T(beta0,x1=x1,x2=x2)/2),pch="+") } # konec pro derivL cat("vzdálenost AB, param. theta", sqrt(crossprod(etaT0-etaT00)),"\n") print(krivost(theta0,x1=x1,x2=x2)) print(krivostT(beta0,x1=x1,x2=x2)) # savePlot("pokus.eps",type="ps") par(opar) } # pomucka pro NLR, priklad fce0 <- function(theta,x1=1,x2=3){ # vektor funkcnich hodnot f(theta) OUT <- NULL for (i in 1:length(theta)) OUT <- rbind(OUT,c(exp(x1*theta[i]),exp(x2*theta[i]))) OUT } fce1 <- function(theta,x1=1,x2=3){ # vektor prvnich derivaci f(theta) OUT <- NULL for (i in 1:length(theta)) OUT <- rbind(OUT,c(x1*exp(x1*theta[i]),x2*exp(x2*theta[i]))) OUT } fce2 <- function(theta,x1=1,x2=3){ # vektor druhych derivaci f(theta) OUT <- NULL for (i in 1:length(theta)) OUT <- rbind(OUT,c(x1^2*exp(x1*theta[i]),x2^2*exp(x2*theta[i]))) OUT } fce2T <- function(theta,x1=1,x2=3){ # průmět druhých derivací do nadroviny # predpoklada skalarni theta OUT <- c(crossprod(c(fce1(theta,x1=x1,x2=x2)),c(fce2(theta,x1=x1,x2=x2))))/ c(crossprod(c(fce1(theta,x1=x1,x2=x2))))*fce1(theta,x1=x1,x2=x2) OUT } fce2N <- function(theta,x1=1,x2=3){ # průmět kolmý na tecnou nadrovinu # predpoklada sklarni theta OUT <- fce2(theta,x1=x1,x2=x2)-fce2T(theta,x1=x1,x2=x2) OUT } krivost <- function(theta,x1=1,x2=3){ # krivosti v parametrizaci theta K.N <- sqrt(crossprod(c(fce2N(theta,x1=x1,x2=x2))))/ crossprod(c(fce1(theta,x1=x1,x2=x2))) K.T <- sqrt(crossprod(c(fce2T(theta,x1=x1,x2=x2))))/ crossprod(c(fce1(theta,x1=x1,x2=x2))) OUT <- c(K.N=K.N,K.T=K.T) OUT } fceT0 <- function(beta,x1=1,x2=3){ # vektor funkcnich hodnot f(beta) OUT <- NULL for (i in 1:length(beta)) OUT <- rbind(OUT,c(beta[i]^x1,beta[i]^x2)) OUT } fceT1 <- function(beta,x1=1,x2=3){ # vektor prvnich derivaci f(beta) OUT <- NULL for (i in 1:length(beta)) OUT <- rbind(OUT,c(x1*beta[i]^(x1-1),x2*beta[i]^(x2-1))) OUT } fceT2 <- function(beta,x1=1,x2=3){ # vektor druhych derivaci f(beta) OUT <- NULL for (i in 1:length(beta)) OUT <- rbind(OUT,c(x1*(x1-1)*beta[i]^(x1-2),x2*(x2-1)*beta[i]^(x2-2))) OUT } fceT2T <- function(beta,x1=1,x2=3){ # průmět druhých derivací do nadroviny # predpoklada skalarni theta OUT <- c(crossprod(c(fceT1(beta,x1=x1,x2=x2)),c(fceT2(beta,x1=x1,x2=x2))))/ c(crossprod(c(fceT1(beta,x1=x1,x2=x2))))*fceT1(beta,x1=x1,x2=x2) OUT } fceT2N <- function(beta,x1=1,x2=3){ # průmět kolmý na tecnou nadrovinu # predpoklada sklarni theta OUT <- fceT2(beta,x1=x1,x2=x2)-fceT2T(beta,x1=x1,x2=x2) OUT } krivostT <- function(beta,x1=1,x2=3){ # krivosti v parametrizaci theta K.N <- sqrt(crossprod(c(fceT2N(beta,x1=x1,x2=x2))))/ crossprod(c(fceT1(beta,x1=x1,x2=x2))) K.T <- sqrt(crossprod(c(fceT2T(beta,x1=x1,x2=x2))))/ crossprod(c(fceT1(beta,x1=x1,x2=x2))) OUT <- c(K.N=K.N,K.T=K.T) OUT }