0

I want to use nls to fit some experimental data to this equation:

Model <- as.formula (F ~  (dhg*K1*G+dhg2*K1*K2*G^2)/(1 + K1*G+K1*K2*G^2)  - d0)

but G has been obtained from:

K1*K2*G^3*+K1*(2*K2*Ht-K2*Gt.M+1)*G^2+K1*(Ht-Gt.M)+1*G-Gt = 0

Is there a way to include G as a function in the formula (so it can be solved in each iteration, for example, by getting the root with polyroot)? For example

G <- function(K1, K2, Ht, Gt){
    D <- K1*K2
    C <- K1*(2*K2*Ht-K2*Gt+1)
    B <- K1*(Ht-Gt)+1
    A <- -Gt
    G <- Re(polyroot(z=c(A,B,C,D)))[which(Re(polyroot(z=c(A,B,C,D))) > 0)]
    return(G)
}

I didn't manage to include G in the formula

Edit

Thanks to Onyambu answer, I made the following working example:

Data <- data.frame(
    Gt.M =c(0.000000000,0.001085104,0.002170209,0.003255313,0.004340418,0.005425522,0.006510627,0.007595731,0.008680836,0.009765940,0.010851045,0.011936149,0.013021254,0.014106358,0.015191463,0.016276567),
    Fitted_d.max =c(7.167307, 7.162498, 7.158839, 7.155603, 7.152150, 7.149008, 7.145748, 7.142980,7.139988, 7.137311, 7.134371, 7.131921, 7.129423, 7.127004, 7.124224, 7.121391))

G <- function(K1, K2, Ht, Gt.M){
    D <- K1*K2
    C <- K1*(2*K2*Ht-K2*Gt.M+1)
    B <- K1*(Ht-Gt.M)+1
    A <- -Gt.M
    G <- Re(polyroot(z=c(A,B,C,D)))[which(Re(polyroot(z=c(A,B,C,D))) > 0)]
    return(G)
}

Model2_1 <- as.formula (Fitted_d.max ~ d0 + (dhg*K1*G(K1, K2, Gt.M,Ht=1E-3)+dhg2*K1*K2*G(K1, K2, Gt.M,Ht=1E-3)^2)/(1 + K1*G(K1, K2, Gt.M,Ht=1E-3)+K1*K2*G(K1, K2, Gt.M,Ht=1E-3)^2))

nls(Model2_1,
        data = Data, 
        start = list(K1 = 1, K2=1, dhg=0.1, dhg2=0.1,d0=7.3))

However it returns errors

Error in nls(Model2_1, data = Data, start = list(K1 = 1, K2 = 1, dhg = 0.1,  : 
  step factor 0.000488281 reduced below 'minFactor' of 0.000976562
In addition: There were 50 or more warnings (use warnings() to see the first 50)

and using the preview() function in package nlstools, it returns:

Error in data.frame(..., check.names = FALSE) : 
  arguments imply differing number of rows: 16, 23

I'm afraid I'm not sure where the problem is

thanks!

edit2

I got it working by changing how G is calculated

library('VGAMextra')

Ft <- function(K1, K2, Gt.M,Ht,d0,dhg,dhg2){
G <- function(K1, K2, Gt.M,Ht){
F <- function(K1, K2, Gt.M,G,Ht) K1*K2*G^3 + K1*(2*K2*Ht-K2*Gt.M+1)*G^2 + (K1*(Ht-Gt.M)+1)* G -Gt.M
Fdev <- function(K1, K2, Gt.M,G,Ht) 3*K1*K2*G^2 + 2*K1*(2*K2*Ht-K2*Gt.M+1)*G + (K1*(Ht-Gt.M)+1)
newtonRaphson.basic(f = F, fprime  = Fdev, a = 0, b= Gt.M, K1=K1, K2=K2, Gt.M=Gt.M,Ht=Ht,nmax=100)}
Gc <- G(K1=K1, K2=K2, Gt.M = Gt.M, Ht=Ht)
Fitted_d.max <- d0 + (dhg*K1*Gc+dhg2*K1*K2*Gc^2)/(1 + K1*Gc+K1*K2*Gc^2)
}

 nls(Fitted_d.max~Ft(K1, K2, Ht=0.5E-3, Gt.M, d0 = Data[1,2],dhg,dhg2),
        data = Data, 
        start = list(K1 = 1500, K2=5, dhg=-0.1, dhg2=-0.19),
        lower= c(0,0, -1,-1),
        trace=TRUE, algorithm="port")

seems to be doing the thing.

Alejandro Méndez
  • 1
  • 2

0 Answers0