## Fitting a numerical ODE solution to data

| categories: data analysis | View Comments

fit_ode

## Fitting a numerical ODE solution to data

function main

clc; clear all; close all


Suppose we know the concentration of A follows this differential equation: , and we have the following data:

Lets make these column vectors

tspan = [0 0.1 0.2 0.4 0.8 1]';
Ca_data = [2.0081  1.5512  1.1903  0.7160  0.2562  0.1495]';
plot(tspan,Ca_data,'ko ') In this case, one could solve the simple ODE, and transform the data to get a simple linear fit. Instead, we are going to develop a method to directly fit the ODE solution to the data.

parguess = [1.3]; % nonlinear fits need initial guesses
[pars, resid, J] = nlinfit(tspan,Ca_data,@myfunc,parguess)

pars =

2.5889

resid =

0
0.0011
-0.0062
0.0031
0.0031
-0.0013

J =

0
-0.1550
-0.2393
-0.2852
-0.2025
-0.1508



the residuals are pretty low. J is a Jacobian, and we don't need to worry about it. It is used to estimate confidence intervals.

## parameter confidence intervals

alpha = 0.05; % this is for 95% confidence intervals
pars_ci = nlparci(pars,resid,'jacobian',J,'alpha',alpha)

sprintf('The true value of k is between %1.2f and %1.2f at the 95%% confidence level',pars_ci)

pars_ci =

2.5701    2.6077

ans =

The true value of k is between 2.57 and 2.61 at the 95% confidence level



## now we can examine the fit

tfit = linspace(0,1);
fit = myfunc(pars,tfit);
hold all
plot(tfit,fit) You can see the fit looks ok. To summarize, here is a way to directly fit an ODE solution to data using the nonlinear fitting algorithms of Matlab. The advantage of this over a transformation, e.g. by computing the derivative or by plotting log(Ca/Cao) vs. something is that this approach avoids the nonlinear transformation of errors in the data. That comes at the expense of having to solve a nonlinear problem which requires an initial guess. Sometimes, that can be hard.

## function to fit

This is the function we want to fit to the data. When we have the right pars, the output of this function should be a close fit to the data above. We use a nested function in a nested function to avoid needing lots of function handles to pass parameters around.

    function output = myfunc(pars,tspan)
% we just rename the variable here for use in the ODE
k = pars;
% this is the ODE we are fitting to

end % main