Re: nonlinear programming with differential-algebraic constraints (2)

*To*: mathgroup at smc.vnet.net*Subject*: [mg52129] Re: [mg52084] nonlinear programming with differential-algebraic constraints (2)*From*: DrBob <drbob at bigfoot.com>*Date*: Thu, 11 Nov 2004 04:53:09 -0500 (EST)*References*: <200411100945.EAA11228@smc.vnet.net>*Reply-to*: drbob at bigfoot.com*Sender*: owner-wri-mathgroup at wolfram.com

>> The main problem seems to be that the differential >> equation does not get evaluated in each step but only once. model (as written) has the same value every time, so why do you want it evaluated more than once? I think you wanted it to depend on "a", but you've written it with a->0.05 hardwired in. I think this does what you want: f[t_] := Exp[-0.1*t] + Random[Real, {0, 0.1}]; data = {#, f@#} & /@ Table[t, {t, 0, 10, 1}]; {xData, yData} = Transpose@data; Clear[model, square] model[(a_)?NumericQ] := y /. First[NDSolve[ {Derivative[1][y][t] == (-a)*y[t], y[0] == 1}, y, {t, 0, 10}]] square[a_] := (#1 . #1 & )[ yData - model[a] /@ xData] NMinimize[{square[a], a >= 0}, a] {0.00638361, {a -> 0.0818508}} Here's a more direct (but fully equivalent) approach: Clear[f, model, square] model[a_][t_] := Exp[(-a)*t] f[t_] := model[0.1][t] + Random[Real, {0, 0.1}] data = Table[{t, f[t]}, {t, 0, 10, 1}]; {xData, yData} = Transpose[data]; square[a_] := (#1 . #1 & )[ yData - model[a] /@ xData] NMinimize[{square[a], a >= 0}, a] This always underestimates a, however, because the noise you've added is always positive, and the error function doesn't take that into account. You can do better this way: Clear[f, model, square] model[a_][t_] := Exp[-a*t] f[t_] := model[0.1]@t + 0.1Random[] {xData, yData} = Transpose@Table[{t, f@t}, {t, 0, 10, 1}]; square[a_] := #.# &[yData - a/2 - model[a] /@ xData] NMinimize[{square@a, a >= 0}, a] {0.0139239, {a -> 0.0999625}} and this does pretty well with no minimization at all: Clear[f, model, inverse] model[a_, t_] = Exp[-a*t]; inverse[a_]@{x_, y_} := -Log[y - .5a]/x f[t_] := model[0.1, t] + 0.1Random[] data = Table[{t, f@t}, {t, 1, 10}]; Mean[inverse[0.1] /@ data] 0.100146 I didn't put a lot of "science" into that. I'm just subtracting the average bias and solving for a in y - .5 a == Exp[-a x]. The LHS of that equation is an unbiased estimate of y (given a and x), but the result isn't necessarily an unbiased estimate of a (given y and x). Bobby On Wed, 10 Nov 2004 04:45:25 -0500 (EST), Joerg Schaber <schaber at molgen.mpg.de> wrote: > I made a little example to illustrate the problem: > > I want to fit a model which is stated as a system of differenttial > equations to some data. > > (* curve to be fitted *) > f[t_]:=Exp[-0.1*t]+Random[Real,{0,0.1}]; > data={#,Evaluate[f[#]]}&/@Table[t,{t,0,10,1}]; > > ListPlot[data,PlotStyle\[Rule]{RGBColor[0,0.2,0.8],PointSize[0.02]}]; > > (* the parameter 0.1 is to be found *) > (* the model as a differential equations *) > model:=NDSolve[{y'[t]==-a*y[t],y[0]==1}/.a->0.05,y,{t,0,10}]; > > (* minimizing sum of squares *) > NMinimize[{squares=Plus@@Table[(Table[{#,y[t]/.Evaluate[model][[1]]/.t->#}&/@ > Table[t,{t, 0, 10, 1}]][[i, 2]]-data[[i,2]])^2, > {i, Length[data]}], a > 0}, a, > Method -> {"DifferentialEvolution"}, > StepMonitor :> Print[squares, a], > WorkingPrecision -> 8] > > > The main problem seems to be that the differential equation does not get > evaluated in each step but only once. > > So if anybody could make this example work that would be great. > > best, > > joerg > > > > -- DrBob at bigfoot.com www.eclecticdreams.net

**References**:**nonlinear programming with differential-algebraic constraints (2)***From:*Joerg Schaber <schaber@molgen.mpg.de>