Fitting of coupled differential equations

Share icon

Hello,

Several places in the Igor manual refer to the fitting of coupled, first order ODEs, and I'm trying to do that. This particular system has been solved analytically, but I imagine future systems where I will not be able to use an analytical solution.

Ultimately, it would be great if this could run from the curve-fitting dialog, but as I will have many datasets to fit, I'd like to have it easily scripted, as well.

I am using IntegrateODE to perform the integration (over time); initially, I had a separate wave of time points where observations had been made for each of the four observables (Maa, Mab, Mbb, Mba). The observables themselves were stored in a 2D wave. This worked fine, but when trying to write a fitting function (FitFunc) realized two things:

- the curve fitting window would not accept a 2D wave as a dependent variable (not critical, but unfortunate)
- when running from the command line, I got the error: "The data must have one dimension for each independent variable"

After doing this, I modified my ODE function and my FitFunc to accept a "master" matrix of dependent variables, where the desired time points are stored in the matrix itself. I realize this may not be the correct solution. Regardless, the code is below:

This is the ODE function:

Function ZZ(pw, xt, yw, dydx)
    Wave pw         // kinetic rate constants
    Variable xt         // mixing time
    Wave yw         // peak intensities corresponding to Maa, Mab, Mbb, and Mba
    Wave dydx       // derivative list (dMaa/dt, dMab/dt, ...)

// parameter indices are (0-1 not used in this function):
//  
//  0   Maa_0       // Initial height of Maa
//  1   Mbb_0       // Initial height of Mbb
//  2   Raa/ba      // Relaxation rate of a and b states
//  3   Rbb/ab      //

//  4   kab     // The (asymmetric) rates -- what you probably actually want
//  5   kba     //
   
    dydx[0] = 0
    dydx[1] = -(pw[2]+pw[4])    * yw[1] +       pw[5]   * yw[2]
    dydx[2] =   pw[4]       * yw[1] +   -(pw[3]+pw[5])  * yw[2]
    dydx[3] = -(pw[3]+pw[5])    * yw[3] +       pw[4]   * yw[4]
    dydx[4] =   pw[5]       * yw[3] +   -(pw[2]+pw[4])  * yw[4]
   
    return 0
End


And here is the fitting function:

Function FitZZ_asym(pw, dataw0, xw) : FitFunc
    Wave pw     // parameter list (values to be fit)
    Wave dataw0     // N by 5 matrix containing mixing times and peak intensities
    Wave xw     // ignored

// parameter indices are desciribed above

// the peak intensities is assumed to be a table with five columns,
// in order the columns correspond to t, Maa, Mab, Mbb, Mba

// Obviously, including the zero point explicitly in the heights will be much
// faster; however, most will not have implemented their table in this way.

    // zero-point included explicitly
    if (dataw0[0][0] == 0)
        dataw0[0][1] = pw[0]
        dataw0[0][2] = 0.0
        dataw0[0][3] = pw[1]
        dataw0[0][4] = 0.0
 
        // IntegrateODE will overwrite the time values in dataw0; this
        // saves them for later    
        Make/O/D/N=(dimsize(dataw0, 0)) tw
        tw = dataw0[p][0]
       
        IntegrateODE /M=1 /X=tw ZZ, pw, dataw0
       
        // put back the time values
        dataw0[][0] = tw
        return 0
    endif
   
    // zero-point not included; must copy and resize
    Duplicate/O/D dataw0 dataw
    InsertPoints 0, 1, dataw
   
    dataw[0][0] = 0.0
    dataw[0][1]     = pw[0]
    dataw[0][2]     = 0.0
    dataw[0][3]     = pw[1]
    dataw[0][4]     = 0.0
   
    // IntegrateODE will overwrite the time values in dataw0; this
    // saves them for later    
    Make/O/D/N=(dimsize(dataw, 0)) tw
    tw = dataw[p][0]
   
    IntegrateODE /M=1 /X=tw ZZ, pw, dataw      

    // put back the time values
    dataw[][0] = tw
   
    // resize new matrix and copy the data back to the original
    DeletePoints 0, 1, dataw
    dataw0 = dataw
   
    return 0
End


Now, if I make a new wave containing some sample data (generated elsewhere) and initial parameters:

Make/D/O/N=(14, 5)
hgt2[0][0]= {0.005,0.01,0.015,0.02,0.05,0.1,0.15,0.2,0.3,0.5,0.75,1,2,5}
hgt2[0][1]= {0.992,0.98,0.97,0.964,0.905,0.825,0.75,0.682,0.574,0.422,0.305,0.225,0.101,0.012}
hgt2[0][2]= {0.005,0.013,0.019,0.029,0.071,0.13,0.172,0.211,0.27,0.329,0.335,0.311,0.175,0.021}
hgt2[0][3]= {1.289,1.281,1.272,1.261,1.207,1.128,1.051,0.983,0.866,0.691,0.54,0.431,0.2,0.027}
hgt2[0][4]= {0.004,0.01,0.014,0.018,0.043,0.079,0.104,0.127,0.164,0.198,0.205,0.187,0.104,0.015}

Make/D/O/N=6 pw0 = {1,1,3,3,0.4,0.4}


I hope to run my fitting routing with a command similar to the following:

FuncFit FitZZ_asym, pw0, hgt2


But here again, I'm getting the "The data must have one dimension for each independent variable" error I got before.

Clearly, I'm not understanding something about how multivariable fits are done. Ultimately, I have to give this code to graduate students so I'd prefer to have the minimum amount of "kludge" so they can understand what's going on. I'd appreciate any advice on what I'm doing wrong here.

Thanks,
Nick
Hi Nicholas,

nfitzkee wrote:


- the curve fitting window would not accept a 2D wave as a dependent variable (not critical, but unfortunate)



If you want  IntegrateODE to produce the integrated values at specific time values you need to provide the time wave (tw) a matrix with a column for each differential equation (sim_hgt) and their respective initial values. The wave tw and the matrix must also have the same number of columns

I would change the differential equation to:

Function ZZ(pw, xt, yw, dydx)
    Wave pw         // kinetic rate constants
    Variable xt         // mixing time
    Wave yw         // peak intensities corresponding to Maa, Mab, Mbb, and Mba
    Wave dydx       // derivative list (dMaa/dt, dMab/dt, ...)
 
// parameter indices are (0-1 not used in this function):
//  
//  0   Maa_0       // Initial height of Maa
//  1   Mbb_0       // Initial height of Mbb
//  2   Raa/ba      // Relaxation rate of a and b states
//  3   Rbb/ab      //
 
//  4   kab     // The (asymmetric) rates -- what you probably actually want
//  5   kba     //
 

    dydx[0] = -(pw[2]+pw[4])    * yw[0] +       pw[5]       * yw[1]
    dydx[1] =   pw[4]       * yw[0] +   -(pw[3]+pw[5])  * yw[1]
    dydx[2] = -(pw[3]+pw[5])    * yw[2] +       pw[4]       * yw[3]
    dydx[3] =   pw[5]       * yw[2] +   -(pw[2]+pw[4])  * yw[3]
 
    return 0
End


and I would run the following commands to create the simulated values.

Make/D/O/N=6 pw0 = {1,1,3,3,0.4,0.4}
Make/D/O/N=15 tw= {0, 0.005,0.01,0.015,0.02,0.05,0.1,0.15,0.2,0.3,0.5,0.75,1,2,5}
Make/D/O/N=(15, 4) sim_hgt
sim_hgt[][0] = {{pw0[0]},{0},{pw0[1]},{0}}
IntegrateODE /M=1 /X=tw ZZ, pw0, sim_hgt


nfitzkee wrote:


The data must have one dimension for each independent variable



This error happens because you are trying to fit four data sets simultaneously. The only "easy" way I know of doing this is to perform a global fit for the data with a fitting function for each column.

For example, in order to fit the first column (i.e. Maa), you would need a function like this.

Function FitZZ_Maa(pw, dataw, xw) : FitFunc
    Wave pw     // parameter list (values to be fit)
    Wave dataw  // N by 4 matrix containing mixing times and peak intensities
    Wave xw     // ignored
 
        Make/O/D/N=(dimsize(dataw, 0),4) sim_data

        sim_data[0][0] = pw[0]
        sim_data[0][1] = 0.0
        sim_data[0][2] = pw[1]
        sim_data[0][3] = 0.0
 
        IntegrateODE /M=1 /X=xw ZZ, pw, sim_data
        dataw = sim_data[p][0]
 
    return 0
End

and you would need a command like this to run the fitting function
FuncFit FitZZ_Maa, pw0, hgt2[][0] /X=tw


nfitzkee wrote:

Ultimately, I have to give this code to graduate students so I'd prefer to have the minimum amount of "kludge" so they can understand what's going on.


I am a graduate student myself :), so I hope this points you in the right direction.

March 25, 2013 at 09:00 am - Permalink

The solution of using Global Fit for this problem is convenient, but will involve solving the ODE system four times for every time the fit needs the solution. It might be possible to figure out a way to cache the first result so it can be looked up later by the other fit functions, but it could be tricky to get it to work correctly.

I talked to the original poster on the phone shortly after he posted- I suggested a scheme in which the four data sets are simply concatenated. This requires that the fitting function be aware of this, and do the right thing with each of the four columns of the solution.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

March 25, 2013 at 10:31 am - Permalink

johnweeks wrote:
The solution of using Global Fit for this problem is convenient, but will involve solving the ODE system four times for every time the fit needs the solution. It might be possible to figure out a way to cache the first result so it can be looked up later by the other fit functions, but it could be tricky to get it to work correctly.

I talked to the original poster on the phone shortly after he posted- I suggested a scheme in which the four data sets are simply concatenated. This requires that the fitting function be aware of this, and do the right thing with each of the four columns of the solution.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com


Exactly how? I was doing a similar fit where a coupled differential equation to be solved and then fit it to two different curve. What I did, I added two waves to form a single wave and did the same for the x waves (its just repeating of the same x wave twice in the new x wave). Then I was trying via all at once type fitting. But its not working in several levels. Can you please help?

July 8, 2016 at 07:21 am - Permalink

It was a long time ago and I don't have any memory of the phone call. But your brief description sounds like the right approach. Naturally, from your information I can't tell what's going wrong or whether you have understood the requirements of the method correctly.

You could contact us at support@wavemetrics.com.

John Weeks
WaveMetrics, Inc.
support@wavemetrics.com

July 8, 2016 at 11:05 am - Permalink