Parameter fitting for a set of rate equations
chriswo
my name is Chris, I am a novice in Igor programming, though I try to improve my skills.
Currently, I am trying to teach Igor how to fit to a set of data (t,y- points), that are somehow related to a set of differential equations, in order to extract fitting parameters.
For learning purposes I checked the simple case with one differential equation, but I am not exactly sure, what is happening, at least not at all stages. (see: http://www.igorexchange.com/node/1940)
I copy the final code and try to comment how I understand it and what the purpose is, please correct me if I am mistaking and help me understand!
Function ExpDecay (pw,xx,yw,dydx)// Definition of a function with dependencies pw, xx, dydx, yw ?
Wave pw, yw, dydx //Here he says: pw, yw and dydx are actual waves and subseq. xx is the variable
Variable xx
dydx[0] = -pw[0]*yw[0] //this is the differential equation, though I am not exactly sure, why one needs the [0] especially the pw[0]. I understand this as the "rate" (as it's called later, but there it becomes the y offset???)
End
Function FitExpDecay(pw, yw, xw) : FitFunc // Here the actual fit-function is defined with parameters (pw) and data points yw,xw
wave pw, yw, xw
//pw[0]= Y shift (y0)
//pw[1]= Scalar (Ao)
//pw[2]= Rate (a)
Make/O/D/n=1 pODE // Ok, what happens here? A wave is created with one column, what for?
pODE[0]=pw[2] // pODE[0] (<-- first entry??) is defined to be equal to pw[2](<-- 3rd entry in the parameter wave??) why?
yw[0]=pw[1] // no clue what this means
IntegrateODE /X=xw ExpDecay, pODE, yw // ok this is the command for the fitting, right? i thought here were only 3 "variables" : ItegrateODE derivFunc(eExpDecay), cwaveName (parameters??), output
yw += pw[0] // no clue, it must somehow relate to a y-offset, since this doesn't appear anywhere else
End
Wave pw, yw, dydx //Here he says: pw, yw and dydx are actual waves and subseq. xx is the variable
Variable xx
dydx[0] = -pw[0]*yw[0] //this is the differential equation, though I am not exactly sure, why one needs the [0] especially the pw[0]. I understand this as the "rate" (as it's called later, but there it becomes the y offset???)
End
Function FitExpDecay(pw, yw, xw) : FitFunc // Here the actual fit-function is defined with parameters (pw) and data points yw,xw
wave pw, yw, xw
//pw[0]= Y shift (y0)
//pw[1]= Scalar (Ao)
//pw[2]= Rate (a)
Make/O/D/n=1 pODE // Ok, what happens here? A wave is created with one column, what for?
pODE[0]=pw[2] // pODE[0] (<-- first entry??) is defined to be equal to pw[2](<-- 3rd entry in the parameter wave??) why?
yw[0]=pw[1] // no clue what this means
IntegrateODE /X=xw ExpDecay, pODE, yw // ok this is the command for the fitting, right? i thought here were only 3 "variables" : ItegrateODE derivFunc(eExpDecay), cwaveName (parameters??), output
yw += pw[0] // no clue, it must somehow relate to a y-offset, since this doesn't appear anywhere else
End
What I actually want to do is a different system/ set of rate equations: something like a two-level system (with no analytic solution) in quantum mechanics:
again I will have a set of (y,t) data points.
the differential equations are then something like:
dN2/dt = N1*k1-N2*(k2+k3)
dN1/dt =-dN2/dt
and they are connected to the data points with:
y(t)=N1(t)+N2(t)
I tried to start like this:
Function NonExpDecay (pw2,xx,yw,dN2dt, dN1dt,N1,N2)
Wave pw2, yw, dN2dt, dN1dt,N1,N2
Variable xx
dN2dt[0] = N1[0]*pw2[0]-N2[0]*(pw2[1]+pw2[2])
dN1dt[0] =-dN2dt[0]
yw[0]=N1[0]*N2[0]
End
Wave pw2, yw, dN2dt, dN1dt,N1,N2
Variable xx
dN2dt[0] = N1[0]*pw2[0]-N2[0]*(pw2[1]+pw2[2])
dN1dt[0] =-dN2dt[0]
yw[0]=N1[0]*N2[0]
End
Is this correct? Again I am not eactly sure about the bookkeeping (the [0]'s ?)
Now I don't know how to proceed?
Of course I will need to introduce a fit-function, something like
Function FitRates(pw2, yw, tw) : FitFunc
wave pw2, yw, tw
//pw[0]= k1 : transition rate Absorption
//pw[1]= k2 : transition rate rad. emission
//pw[2]= k3 : transition rate non-rad. emission
.....
wave pw2, yw, tw
//pw[0]= k1 : transition rate Absorption
//pw[1]= k2 : transition rate rad. emission
//pw[2]= k3 : transition rate non-rad. emission
.....
Can somebody please help?
Thank you in advance!!
Best wishes,
Chris
Have you read through the help on Solving Differential Equations? Copy this command, paste it into Igor's command line and press Enter:
DisplayHelpTopic "Solving Differential Equations"
The simple example you show is included in that help file in the sub-topic "A First-Order Equation". That help file goes on to examples showing a system of first-order equations (that's what you want to solve) and then on to a second-order example.
Big Tip: Take this one step at a time. Understand the help on the first-order equation, and the system of first-order equations. When you can make the examples work, then try your own equation. Only at that point should you even think about trying to incorporate the differential equation solution into a fit function.
So, read your homework. Post further questions :)
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
August 7, 2014 at 02:25 pm - Permalink
wave pw, yw, xw
//pw[0]= Y shift (y0)
//pw[1]= Scalar (Ao)
//pw[2]= Rate (a)
Make/O/D/n=1 pODE // Generate a wave for one of the output parameters
pODE[0]=pw[2] // here it is actually defined: the rate a
yw[0]=pw[1] // Ok so here one links the fitting values into the fitting function. The first, initial value is put in.
IntegrateODE /X=xw ExpDecay, pODE, yw // so now, this command says: integrate with xw as wave using the function ExpDecay, give back the parameters in pODE and the solution in yw
yw += pw[0] // isn't a certain knowledge about the solution necessary, though I understand, that a y-offset is not really specific, but still?
End
Back to my example:
I guess one can make this a bit nicer:
Wave pw, yw, dydt
Variable tw
dydt[0] = yw[0]*pw[0]-yw[1]*(pw[1]+pw[2])
dydt[1] =-dydt[0]
yw[2]=yw[0]*yw[1] // is this going to work? One could certainly define this as a differential equation: dydt[2]=dydt[0]*yw[1]+dydt[1]*yw[0]
End
Now I want to define the fit function such, that it fits to yw[2], so in principle the product of the solutions of the two differential equations.
And in particular this gives me a hard time.
Wave pw // parameters
Wave dataw // N by 3 matrix
Wave xw // ignored
Make/O/D/N=(dimsize(dataw, 0),3) fitdata
fitdata[0][0] = pw[0]
fitdata[0][1] = pw[1]
fitdata[0][2] = pw[2]
IntegrateODE /M=1 /X=xw twolevel, pw, fitdata
dataw = fitdata[p][2] //in order to specify which data points I actually want to fit?
return 0
End
Oh yeah and certainly:
FuncFit fitPL, para, myactualYdata /X=myactualXdata
Is this going to work?
Thanks in advance!!
Best wishes, Chris
September 16, 2014 at 01:53 am - Permalink
You're getting there.
What you have is a second-order system involving two state variables (N1 and N2), plus an equation that connects those two state variables to your observed variable. IntegrateODE will solve for N1 and N2; you will have to evaluate that last equation separately as a post-processing step after the ODE system is solved.
So the output of IntegrateODE will be a two-column wave, with one column for each of N1 and N2. Then after IntegrateODE is done, you will add the two columns together.
So remove the last line from twolevel(); that needs to be done later. Put it into the fitting function. The derivative function, twolevel(), will be called many, many times to evaluate the derivatives at a bunch of points in the space of the ODE system. There is nothing in twolevel() that actually knows about your final state. The calls aren't in order in any sense, and aren't necessarily on the actual final trajectory of the solution. The solver has to make some calls just to find out enough information about the system to be able to tiptoe through the space looking for a trajectory that satisfies the error conditions and solves the equations.
Please don't try doing a curve fit until you have successfully run IntegrateODE... the fitting function can serve as a way to drive IntegrateODE, but actually doing the fit will put an additional level of complexity on top of figuring out how to solve the ODEs.
In your FitPL() function, you almost have the idea. But your data, if I understand correctly, will be 1D- that is, one column. It is to be fit to a combination of the state variables in your ODE system. That, in turn, means that you need to create an intermediate wave with two columns (right? I'm not sure why you think you need three columns in FitPL()) and then post-process the two columns into the final form to fit your data. Something like this:
Wave pw // parameters
Wave dataw // Just N
Wave xw // The actual X values where you need solution points; not ignored!
Make/FREE/D/N=(numpnts(dataw)) ODESolution
ODESolution[0][0] = ??? initial condition for N1, maybe involving elements of pw if the initial conditions are to be fitted
ODESolution[0][1] = ??? initial condition for N2, maybe involving elements of pw if the initial conditions are to be fitted
IntegrateODE /M=1 /X=xw twolevel, pw, ODESolution
dataw = ODESolution[p][0] + ODESolution[p][1]
return 0
End
If you don't understand what's happening in the assignment to dataw, please read the help:
DisplayHelpTopic "Waveform Arithmetic and Assignment"
DisplayHelpTopic "Multidimensional Wave Indexing"
DisplayHelpTopic "Multidimensional Wave Assignment"
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
September 16, 2014 at 09:38 am - Permalink