Fitting a cubic-like curve to data in R

Question

I tried using nls() in R to fit the following expression to a set of data:

where g=9.8, alpha & B_0 are unknown, a = 0.01, z_0 = 0.3 such that:

theory <- as.formula(V~-(9.8)ab*(pi)(0.02^2)(T)(0.3-(0.5)(9.8)(T^2))exp(-a(0.3-(0.5)(9.8)*(T^2))^2))

where T = t and V = E(t) from the theory model

I was expecting a cubic-like curve but it resulted in a relatively straight line. Is there another way to try and fit the curve, such that the result is closer to:

EDIT: here are the results I got when using nls:

theory <- as.formula(V~-(9.8)ab*(pi)(0.02^2)(T)(0.3-(0.5)(9.8)(T^2))exp(-a(0.3-(0.5)(9.8)*(T^2))^2))

mydata <- read.csv("voltage1.csv") nls.1 <- nls(theory,start=list(a=.8,b=.01),data=mydata) overview(nls.1) plotfit(nls.1,smooth=TRUE)

Answer 1

The two peaks arise from this term $(z-0.5gt^2)$ and its square $(z-0.5gt^2)^2$.

Moving away from the point where $(z-0.5gt^2) = 0$, the peaks increase from the increase in the magnitude of $(z-0.5gt^2)$ (which can be negative/positive) and the peaks decrease when the terms $e^{-(z-0.5gt^2)^2}$ start to become dominant.

In the plot below you see that this term $(z-0.5gt^2)^2$ is zero around the point $t = \sqrt{2z/g} \approx 0.247$. The negative peak will be on the left of it and the positive peak will be on the right of it.

Your data has the both peaks on the left of $t = 0.25$, and that is why there is no good fit. Probably your value of $z$ is wrong or the time $t$ needs to be shifted.

If you allow flexibility for either $z$ or if you let the parameter $t$ be flexible with a shift, then you are able to get a fit that is not a straight line because now there is are two peaks possible. However, compared to the data the peaks are more broad (the red curve in the plot). It can be fixed by using a slightly different formula (the green curve in the plot)

$$v = -g \alpha b_0 \pi a^2 (t) (z-0.5 g (t)^2)^3 e^{-\alpha(z-0.5g(t)^2)^2}$$

The strange thing about your formula without the power $(z-0.5 g (t)^2)^3$ is that the derivative $dv/dt$ is not zero when this term $(z-0.5 g (t)^2)$ is zero (the point in between the two peaks. But in the plots we see that in between the peaks the slope is zero, or at least there is a strong inflection point.

t = c(seq(0,0.15,0.01), 0.155, seq(0.16,0.25,0.01))
v = c(0,-0.1,-0.2,-0.3,-0.4,-0.8,-1.3,-2.6,-2.4,-1.6,-0.5,-0.2,0,0.4,1,3.6,4.6,3.6,1,0.2,rep(0,7))
plot(t,v)

g = 9.8
a = 0.01
z = 0.3
mod = nls(v ~ -g*alpha*b*pi*a^2*(t+t0)*(z-0.5*g*(t+t0)^2)*exp(-alpha*(z-0.5*g*(t+t0)^2)^2) , start = c(alpha = 200, b = 500, t0 = 0.15))
mod3 = nls(v ~ -g*alpha*b*pi*a^2*(t+t0)*(z-0.5*g*(t+t0)^2)^3*exp(-alpha*(z-0.5*g*(t+t0)^2)^2) , start = c(alpha = 200, b = 50000, t0 = 0.15))

ts = seq(0,0.25,0.0001)
lines(ts, predict(mod, newdata = list(t=ts)), col = 2)
lines(ts, predict(mod3, newdata = list(t=ts)), col = 3)