Re: nonlinear programming with differential-algebraic constraints (2)
- To: mathgroup at smc.vnet.net
- Subject: [mg52120] Re: nonlinear programming with differential-algebraic constraints (2)
- From: Peter Pein <petsie at arcor.de>
- Date: Thu, 11 Nov 2004 04:52:37 -0500 (EST)
- References: <cmsok4$b72$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Joerg Schaber 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 > Hi Jörg, I changed the code to let me play more easily with different ranges and counts of points. I changed f to produce a relative error and finally I used FindMinimum instead of NMinimize because I've got no access to version 5 of Mathematica: In[1]:= (*curve to be fitted*) n = 50; f[t_] := Exp[-0.1*t]*Random[Real, {0.9, 1.1}]; data = {#, Evaluate[f[#]]} & /@ (10 Range[0, n]/n); ListPlot[data, PlotStyle -> {RGBColor[0, 0.2, 0.8], PointSize[0.02]}, Axes -> False, Frame -> True]; (*Plot ommitted*) In[5]:= (*the parameter 0.1 is to be found*) (*the model as a differential equations*) (*evaluate only for numeric a, and increment a counter*) model[a_?NumericQ] := (solvecount++; NDSolve[{y'[t] == -a*y[t], y[0] == 1}, y, {t, Sequence @@ First /@ data[[{1, -1}]]}]); In[6]:= solvecount = 0; FindMinimum[(#1 . #1 & )[ (y[#1] /. Evaluate[model[a]][[1]] & ) /@ First /@ data - Last /@ data ], {a, 0, 1}] solvecount Out[7]= {0.0672762, {a -> 0.101174}} Out[8]= 816 As you can see, model[] has been called 816 times. I hope this helped, Peter -- Peter Pein Berlin