#### See Shi et al. (2025) for details #########################################
# Shi, P., Liu, X., Gielis, J., Beirinckx, B., Niklas, K.J. (2025)
# Comparison of six non-linear equations in describing the 2-D
# profiles of apical meristems. American Journal of Botany (under review).
################################################################################
data(SAMs)
uni.sam <- sort( unique(SAMs$Genus) )
ind <- 2
Data <- SAMs[SAMs$Genus==uni.sam[ind], ]
x0 <- Data$x
y0 <- Data$y
Res1 <- adjdata(x0, y0, ub.np=200, times=1.2, len.pro=1/20)
X <- Res1$x
Y <- Res1$y
dev.new()
plot( X, Y, pch=1, cex.lab=1.5, cex.axis=1.5,
xlab=expression(paste(italic(x), " (unitless)", sep="")),
ylab=expression(paste(italic(y), " (unitless)", sep="")) )
x1 <- X
y1 <- Y
if(TRUE){
# The code results in an upward-opening profile curve with
# all y-values less than zero, facilitating data fitting by
# providing suitable initial values for model parameters.
y1 <- min(Y)-Y
x1 <- X-mean(X)
# To normalize the profile coordinates,
# and make x-coordinates range between -1 and 1
x2 <- x1/max(abs(x1))
y2 <- y1/max(abs(x1))
rm(x1)
rm(y1)
x1 <- x2
y1 <- y2
}
#### Fitting the catenary equation ###############################################
myfun1 <- function(P, x){
a <- P[1]
return( a * cosh(x/a) )
}
x0.ini <- 0
y0.ini <- -min(y1)
theta.ini <- 0
a.ini <- c( 0.75*abs(min(y1)), abs(min(y1)), 1.25*abs(min(y1)) )
ini.val1 <- list(x0.ini, y0.ini, theta.ini, a.ini )
Res1 <- fitAM( myfun1, x=x1, y=y1, ini.val=ini.val1,
control=list(reltol=1e-20, maxit=20000),
par.list=FALSE, stand.fig=FALSE, fig.opt=TRUE, angle=NULL,
unit="unitless", main="Catenary equation" )
###################################################################################
#### Fitting the parabolic equation ###############################################
myfun2 <- function(P, x){
alpha1 <- P[1]
alpha2 <- P[2]
alpha1*x + alpha2*x^2
}
x0.ini <- 0
y0.ini <- -abs(min(y1))
theta.ini <- 0
beta1.ini <- 1e-3
beta2.ini <- c(1, 5, 10, 15)
ini.val2 <- list(x0.ini, y0.ini, theta.ini, beta1.ini, beta2.ini)
Res2 <- fitAM( myfun2, x=x1, y=y1, ini.val=ini.val2,
control=list(factr=1e-12, maxit=20000),
method="L-BFGS-B", lower=c(-Inf, -Inf, -pi/2/4, -Inf, 0),
upper=c(Inf, Inf, pi/2/4, Inf, Inf),
par.list=FALSE, stand.fig=FALSE, fig.opt=TRUE, angle=NULL,
unit="unitless", main="Parabolic equation" )
###################################################################################
#### Fitting the hybrid catenary-parabolic equation ###############################
myfun3 <- function(P, x){
alpha <- P[1]
beta <- P[2]
gamma <- P[3]
alpha * cosh(beta * x) + gamma * x^2 - alpha
}
x0.ini <- 0
y0.ini <- -abs(min(y1))
theta.ini <- 0
alpha.ini <- c(0.1, 0.5, 1, 2.5, 5)
beta.ini <- 1/abs(min(y1))
gamma.ini <- 1
ini.val3 <- list(x0.ini, y0.ini, theta.ini, alpha.ini, beta.ini, gamma.ini)
Res3 <- fitAM( myfun3, x=x1, y=y1, ini.val=ini.val3,
control=list(reltol=1e-30, maxit=50000),
par.list=FALSE, stand.fig=FALSE, fig.opt=TRUE, angle=NULL,
unit="unitless", main="Hybrid catenary-parabolic equation" )
###################################################################################
#### Fitting the performance equation #############################################
myfun4 <- function(P, x){
c <- P[1]
K1 <- P[2]
K2 <- P[3]
xmin <- 0
xmax <- P[4]
# x[x > xmax] <- xmax
# x[x < xmin] <- xmin
-c * ( 1-exp(-K1*(x-xmin)) )*( 1-exp(K2*(x-xmax)) )
}
x0.ini <- min(x1)
y0.ini <- 0
theta.ini <- -pi/12
c.ini <- abs(min(y1))
K1.ini <- 1
K2.ini <- 1
xmax.ini <- max(x1)*2
ini.val4 <- list(x0.ini, y0.ini, theta.ini, c.ini, K1.ini, K2.ini, xmax.ini)
Res4 <- fitAM( myfun4, x=x1, y=y1, ini.val=ini.val4, method="L-BFGS-B",
lower=c(-Inf, -Inf, -pi/2/4, 0, 0, 0, 0),
upper=c(Inf, Inf, pi/2/4, Inf, Inf, Inf, Inf),
control=list(factr=1e-12, maxit=20000),
par.list=FALSE, stand.fig=FALSE, fig.opt=TRUE, angle=NULL,
unit="unitless", main="Performance equation" )
###################################################################################
#### Fitting the superparabolic equation equation #################################
myfun5 <- function(P, x){
beta1 <- P[1]
beta2 <- P[2]
beta1 * abs(x)^beta2
}
x0.ini <- 0
y0.ini <- -max( y1 )
theta.ini <- 0
beta1.ini <- max(x1)/2
beta2.ini <- c(1.5, 2, 2.5)
ini.val5 <- list(x0.ini, y0.ini, theta.ini, beta1.ini, beta2.ini)
Res5 <- fitAM( myfun5, x=x1, y=y1, ini.val=ini.val5,
control=list(reltol=1e-20, maxit=20000),
par.list=FALSE, stand.fig=FALSE, fig.opt=TRUE, angle=NULL,
unit="unitless", main="Superparabolic equation" )
###################################################################################
#### Fitting the superellipse equation equation ###################################
myfun6 <- function(P, x){
A <- P[1]
B <- P[2]
n <- P[3]
-B*(1-abs(x/A)^n)*(1/n)
}
x0.ini <- 0
y0.ini <- 0
theta.ini <- c(0, pi/2)
A.ini <- max(x1)
B.ini <- -min(y1)
n.ini <- c(1, 2, 3)
ini.val6 <- list(x0.ini, y0.ini, theta.ini, A.ini, B.ini, n.ini)
Res6 <- fitAM( myfun6, x=x1, y=y1, ini.val=ini.val6,
control=list(reltol=1e-20, maxit=20000),
par.list=FALSE, stand.fig = FALSE, fig.opt=TRUE, angle=NULL,
unit="unitless", main="Superellipse equation" )
###################################################################################
graphics.off()
Run the code above in your browser using DataLab