[R-sig-dyn-mod] iterative procedure vs deSolve

Johannes Ranke johannes.ranke at jrwb.de
Thu Mar 17 18:26:03 CET 2016 

Dear Paula,

did you try to adapt rtol and atol, which are passed from ode() to its default 
solver lsoda(), to lower values? I had a similar case when I compared an 
Eigenvalue based solution of a system of ODEs with the solution produced by 
deSolve, and it turned out the the deSolve solution converged to the 
Eigenvalue based solution when atol and/or rtol were decreased.

Kind regards,

Johannes


Am Donnerstag, 17. März 2016, 14:56:40 schrieb Paula Sanchez Marin:
> Thank you Daniel,
> 
> The analytical solution of the ODE is a system of two non-linear
> equations with two variables.
> You say I can use those expressions and evaluate them at my time points
> with the parameter values and compare the result with the output from ode()
> That's what I was trying to do in excel. The way I know of solving this
> is by using an iterative procedure. If there is another way I am missing
> I will be glad to know!
> 
> Since Crest is needed to calculate Cgills and Cgills is needed to
> calculate Crest, there is a "circular reference" there, that excel
> solves by performing a repeated recalculation of the worksheet cells
> containing the formulas until a specific numeric condition is met (until
> the change is minimal).
> 
> Cheers,
> 
> Paula
> 
> > ------------------------------
> > 
> > Message: 2
> > Date: Thu, 17 Mar 2016 10:24:54 +0100
> > From: Daniel Kaschek <daniel.kaschek at physik.uni-freiburg.de>
> > To: Special Interest Group for Dynamic Simulation Models in R
> > 
> > 	<r-sig-dynamic-models at r-project.org>
> > 
> > Subject: Re: [R-sig-dyn-mod] iterative procedure vs deSolve
> > Message-ID: <1458206694.5529.0 at physik.uni-freiburg.de>
> > Content-Type: text/plain; charset=utf-8; format=flowed
> > 
> > Dear Paula,
> > 
> > I did not really understand what you mean by "iterative procedure" in
> > Excel. You write that it is based on the analytical solution of the ODE:
> > 
> > Cgills=Cgills_ini*exp(-(k01+k12)*t)+(Cu*k01+Crest*k21*wrest/wgills)/(k10+k
> > 12)*(1-exp(-(k10+k12)*t))
> > Crest=Crest_ini*exp(-k21*t)+Cgills*k12*wgills/wrest/k21*(1-exp(-k21*t))
> > 
> > Have you checked the solution? You could use your expressions above and
> > evaluate them at your time points with the parameter values and compare
> > the result with the output from ode(). If they do not coincide, either
> > your analytical solution is wrong or you misspecified the ODE.
> > 
> > Cheers,
> > Daniel
> > 
> > 
> > On Mi, M?r 16, 2016 at 7:59 , Paula Sanchez Marin
> > 
> > <paulasanchez at uvigo.es> wrote:
> >> Hello,
> >> 
> >> While trying to fit some experimental data to a pharmacokinetic model
> >> I have observed that the results of the model predicted by deSolve
> >> are not the same as those resulting from plotting of the solved
> >> equations using an iterative procedure.
> >> I would like to know why do the two methods give different results
> >> and which one is the correct one.
> >> 
> >> Here it goes my example:
> >> A contaminant (Cu) enter an organism by compartment 1 "gills" and
> >> from there it can be transferred to compartment 2 "rest".
> >> The contaminant can only be eliminated from compartment 1 to the
> >> environment.
> >> The parameters (transfer coefficients) are: k01, k10, k12 and k21.
> >> I introduced the concentration or dilution factors when Cu is
> >> transfered from 1 to 2 and viceversa by using the ratios between the
> >> weight of the gills (wgills) and the weight of the rest (wrest).
> >> 
> >> The differential equations describing changes in concentration in
> >> both compartments with time are:
> >> dCgills/dt = k01*Cu+k21*wrest/wgills*Crest-(k10+k12)*Cgills
> >> dCrest/dt=k12*wgills/wrest*Cgills-k21*Crest
> >> 
> >> This is the script I used for solving the differential equations in R:
> >> 
> >> #Constants
> >> Cu=0.03
> >> wrest=0.18
> >> wgills=0.015
> >> #Model
> >> mod<-function(t,y,parms){
> >> 
> >>    with(as.list(c(parms,y)), {
> >>    
> >>      dCgills<-k01*Cu+k21*wrest/wgills*Crest-(k10+k12)*Cgills
> >>      dCrest<-k12*wgills/wrest*Cgills-k21*Crest
> >>      list(c(dCgills,dCrest))
> >>    
> >>    })
> >> 
> >> }
> >> #Parameters
> >> parms=c(k10=0.211,k12=0.234,k01=1196,k21=0.073)
> >> #Initial conditions
> >> X0<-c(Cgills=5.771,Crest=5.905)
> >> #Output
> >> out<-as.data.frame(lsoda(y=X0,times=0:100,func=mod,parms))
> >> write.csv2(out,"out.csv")
> >> 
> >> I used the same Constants (Cu, wrest and wgills), and the same
> >> Parameters and Initial conditions and calculated in excel (allowing
> >> iterations) the results using the equations below (the solution to
> >> the differential equations):
> >> Cgills=Cgills_ini*exp(-(k01+k12)*t)+(Cu*k01+Crest*k21*wrest/wgills)/(k10+
> >> k12)*(1-exp(-(k10+k12)*t))
> >> Crest=Crest_ini*exp(-k21*t)+Cgills*k12*wgills/wrest/k21*(1-exp(-k21*t))>> 
> >> Here are the results (for the first 10 hours) from R:
> >>      time    Cgills     Crest
> >> 
> >> 1      0   5.77100  5.905000
> >> 2      1  36.80185  5.914667
> >> 3      2  56.88706  6.394786
> >> 4      3  70.21069  7.149393
> >> 5      4  79.34640  8.058456
> >> 6      5  85.87763  9.048616
> >> 7      6  90.77894 10.075197
> >> 8      7  94.65042 11.111150
> >> 9      8  97.86195 12.140248
> >> 10     9 100.64178 13.152906
> >> 11    10 103.13106 14.143612
> >> 
> >> And from the iterative procedure in excel:
> >> t 	Cgills 	Crest
> >> 0 	5,771 	5,905
> >> 1 	37,0318579 	6,18568653
> >> 2 	58,2609877 	7,21694678
> >> 3 	73,4040942 	8,60007128
> >> 4 	84,5860141 	10,1314063
> >> 5 	93,0882873 	11,7033524
> >> 6 	99,7407581 	13,260262
> >> 7 	105,102356 	14,7754342
> >> 8 	109,556482 	16,2380253
> >> 9 	113,367627 	17,6454062
> >> 10 	116,718031 	18,998768
> >> 
> >> 
> >> And here the plot:
> >> 
> >> 
> >> Why do the curves differ?
> >> 
> >> Thank you very much for your guindance.
> >> 
> >> _______________________________________________
> >> R-sig-dynamic-models mailing list
> >> R-sig-dynamic-models at r-project.org
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
> > 
> > ------------------------------
> > 
> > Subject: Digest Footer
> > 
> > _______________________________________________
> > R-sig-dynamic-models mailing list
> > R-sig-dynamic-models at r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models
> > 
> > ------------------------------
> > 
> > End of R-sig-dynamic-models Digest, Vol 101, Issue 1
> > ****************************************************
-- 
PD Dr. Johannes Ranke
Wissenschaftlicher Berater

Kronacher Str. 8
79639 Grenzach-Wyhlen
Germany
http://jrwb.de
USt-IdNr.: DE292883940

More information about the R-sig-dynamic-models mailing list