# Numerical solution of varying coefficients ODE

I have a set of observed raw data and use 2nd order ODE to fit the data

$$y''+b_1(t)y'+b_0(t)y = 0$$

The $b_1(t)$ and $b_0(t)$ are time-dependent and I use principal differential analysis(PDA) (R-package: fda, function: pda.fd)to get the estimate of $\hat{b_1}(t)$ and $\hat{b_0}(t)$ .

To check the validity of the estimates $\hat{b_1}(t)$ and $\hat{b_0}(t)$, I use collocation method (R-package bvpSolve, function:bvpcol) to get the numerical solution of the ODE $$y''+\hat{b_1}(t)y'+\hat{b_0}(t)y = 0$$

Then compare the solution from ODE above with the smoothing curve fitting of the raw data.

My question is that my numerical solution from bvpcol can caputure the shape of the fitting curve but not for the value of the function. There are different in term of some constant multiples.

See the figure of my output. The gray dot is my raw data, the red line is Fourier expansion of the raw data, the green line is numerical solution of bvpcol function and the blue line the green-line/1.62. We can see the green line can capture the shape but with values that are constant times of fourier expansion.

I fit several other data and have similar situation but different constant. I am wondering it is the problem of numerical solution of ODE or some other reasons and how to solve this problem to get a good accordance between numerical solution(green) and true Fourier expansion?

Any help and idea is appreciated!

Here is a raw data and code: RData is here

library(fda)
library(bvpSolve)
tvec = 1:length(y)
tvec = (tvec-min(tvec))/(max(tvec)-min(tvec))
# create basis
fbasis = create.fourier.basis(c(0,1),nbasis=47)
bbasis = create.bspline.basis(c(0,1),norder=8,nbasis=47)
bfdPar = fdPar(bbasis)
yfd = smooth.basis(tvec,y,fbasis)$fd yfdlist = list(yfd) bwtlist = rep(list(bfdPar),2) # PDA fit bwt = pda.fd(yfdlist,bwtlist)$bwtlist
# output of estimated coefficients
beta0.fd<-bwt[[1]]$fd beta1.fd<-bwt[[2]]$fd
# define the vary-coef function in terms of t
fbeta0<-function(t)eval.fd(t,beta0.fd)
fbeta1<-function(t)eval.fd(t,beta1.fd)
# define 2nd order ODE
fun2 <- function(t,y,pars) {
with(as.list(c(y,pars)),{
beta0 = pars[[1]];
beta1 = pars[[2]];
dy1 = y[2]
dy2 = -beta1(t)*y[2]-beta0(t)*y[1]
return(list(c(dy1,dy2)))
})
}
# BVP
yinit<-c(p1[1],NA)
yend<-c(p1[length(p1)],NA)
t<-seq(tvec[1],tvec[length(tvec)],0.005)
col<-bvpcol(yini=yinit,yend=yend,x=t,func=fun2,parms=c(fbeta0,fbeta1),atol=1e-5,islin=T)
# plot output
plot(col[,1],col[,2],col='green',type='l')
points(tvec,p1,col='darkgray')
lines(yfd,col='red',lwd=2)
lines(col[,1],col[,2],col='green',type='l')
lines(col[,1],col[,2]/1.62,col='blue',type='l',lwd=2,lty=4)
legend('topleft',col=c('green','darkgray','red','blue'),
legend=c('ODE solution','raw data','basis curve fitting','ODE solution/1.62'),lty=1)


An updated question: are $\hat{b_1}(t)$ and $\hat{b_0}(t)$ unique?

As I change the basis(say bspline,monomial) for expansion of $b_0(t)$ and $b_1(t)$, usually I have very different estimated curves, but the numerical solution of ODE (green line) are very close. That makes me believe different functional form of $\hat{b_1}(t)$ and $\hat{b_0}(t)$ can contribute the same solution. I started a new post with much more detail here.

This is a linear ODE. If you have Dirichlet BC ($y(0) = y_0$, $y(1)=y_1$) the solution will be of the form $y_0*\text{sol}_0+y_1*\text{sol}_1$. How accurate are your boundary conditions compared to raw data and the basis curve fitting conditions and the boundaries?

Looking closely at the plot, you can see the ODE solution goes through the first and last points in the dataset as it should since those are used as boundary conditions. Whereas the curve fitting solution is lower on the right boundary and higher on the left boundary than the data because it is smoothing an interpolating.

If I change the bvp problem to use a value for y(1) scaled to y(1)/1.62, I get a new curve (solid black line) that overlaps the ODE solution/1.62 curve. (Note I had to take p1=y and nbasis=47 in the included code to run it.)

This is not unexpected since it is a linear ode. The real evidence is from looking at the original plot where the fit goes below the points on the right hand side. Constraining the ode solution to stop at this point is why the whole curve is 1.62 higher. If you use the values of y(0) and y(1) from the basis curve fitting, then your ODE model will match much closer. This also tells you something about the sensitivity of your model.

• While I think this is an attempt to respond to the question and contains several suggestions in response to issues raised in the question, I think this answer could be improved in several ways; in particular, if the OP responds to your questions, it might help if you could expand on your responses. Dec 17 '13 at 2:51
• @user21387 Thanks for your response. As you recommend, if I change the boundary condition to be the basis curve fitting, I got almost exact match. Dec 17 '13 at 20:13
• One more question is the estimation of $b_0(t),b_1(t)$ are unique? Dec 17 '13 at 20:29
• no, the estimates of the bc from the curve fitting will depend on the type of fit (spine, spectral, least squares,etc). You could make a measure of the quality of the fit of the ODE to the data (least squares estimate of RMS error) to characterize the model. It could be there is something about the domain of your problem that would inform your choice or constrain your choice of BC. EG you know at t=0, that y(t=0)=0. Dec 18 '13 at 14:37
• @user21387 Can you elaborate a little bit more on that? I started a new post about this here: stats.stackexchange.com/questions/79969/… Dec 18 '13 at 16:12