MathGroup Archive 1996

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Nonlinear fit to diff eqs.

  • To: mathgroup at smc.vnet.net
  • Subject: [mg5086] Re: Nonlinear fit to diff eqs.
  • From: rubin at msu.edu (Paul A. Rubin)
  • Date: Fri, 25 Oct 1996 22:48:51 -0400
  • Organization: Michigan State University
  • Sender: owner-wri-mathgroup at wolfram.com

In article <53gpue$5am at dragonfly.wolfram.com>,
   mod <c530789 at showme.missouri.edu> wrote:
->Dear All,
->
->My problem is to fit a set of data to a system of coupled first order
->ordinary differential equations. The example is:
->
->   Experimental data of y1, y2 and y3 as a function of time
->             y1 = {a1,a2,a3,a4,.....,aN}
->             y2 = {b1,b2,b3,b4,.....,bN}
->             y3 = {c1,c2,c3,c4,.....,cN}
->             time = {t1,t2,t3,t4,.....,tN}
->             
->    Equations to fit (for example):
->
->             y1'[t] = k1 y1[t] + k2 (1-y1[t]-y2[t]) + k4 y3[t]
->             y2'[t] = k2 y2[t] - k3 (1-y2[t]) + k4 y3[t]
->             y3'[t] = k3 (1-y1[t]-y2[t]-y3[t]) +k2 y3[t] 
->
->I need to fit my data to get those k1, k2, k3, and k4. What is the 
easiest
->to to do this with Mathematica? Any package? Please help.
->
->Thaks in advance,
->
->SF

What I don't know about ODEs fills volumes (and in fact did - it filled my 
textbooks from my undergrad DE courses :-), but since I don't see any 
responses here, I'll try giving you a partial answer.

I assume you want to estimate k1 ... k4 to give the minimum 
sum-of-squared-errors fit to the data.  For starters, we can put the 
equations and functions in lists and define a function to solve the ODEs:

  eqs = {y1'[t] == k1 y1[t] + k2 (1-y1[t]-y2[t]) + k4 y3[t],
         y2'[t] == k2 y2[t] - k3 (1-y2[t]) + k4 y3[t],
         y3'[t] == k3 (1-y1[t]-y2[t]-y3[t]) +k2 y3[t]}; 

  funcs = {y1[t], y2[t], y3[t]};

  f[ k_List ] := 
    Flatten[
      funcs /.
        DSolve[ eqs /. Thread[ {k1, k2, k3, k4} -> N[ k ] ], 
                funcs, t ]
    ] /;
    Length[ k ] == 4

Function f takes a list of four numbers as input, solves the diffeqs, and 
spits up a solution.  I've checked to make sure the list has four entries, 
but I haven't put in a check to make sure they're numbers (easy to add).  
Also, I took the liberty of forcing them to machine precision (N[ k ]) to 
avoid getting some seriously lengthy output as Mma tries to use "infinite 
precision" where possible.

A typical application of f[] would look like this:

  f[ {1, 2, -1, 3} ] /. t -> {-1., 1.}

                                   1
  Power::infy: Infinite expression -- encountered.
                                   0.
<repeated a few times>
  General::stop: 
   Further output of Power::infy
     will be suppressed during this calculation.

  {{2. - 0.00238455 C[1] - 8.26016 C[2] + 0.26013 C[3], 
   2. - 9.83365   C[1] - 0.10936 C[2] + 1.92212 C[3]}, 
  {-0.333333 - 0.010591 C[1] -   1.85976 C[2] - 
    0.26013 C[3], -0.333333 - 43.6763 C[1] - 
      0.0246222 C[2] - 1.92212 C[3]}, 
  {-0.222222 - 0.0111639 C[1] + 1.96036 C[2] - 
              -17
    1.14325 10    C[3], 
   -0.222222 - 46.0389 C[1] + 0.025954 C[2] - 
              -17
    8.44757 10    C[3]}}

  Dimensions[%]
  {3, 2}

In this example, I've fed f[] estimates of the four parameters and 
substituted into the ODE solution two choices of t (-1, +1).  The result is 
a 3x2 matrix in which each row is the vector {y1[t], y2[t], y3[t]} for a 
different value of t.

The basic game plan is to write a function that takes as input the list k 
and spits back the sum of squared differences of observed vs. estimated y 
values.  You would be summing over all time epochs and across the three 
functions y1, y2, y3 to get a single total squared error.  That's the 
function you minimize (with FindMinimum).  That function will use the 
function f[] above to get teh "estimated" y values.  You'll need to make 
those constants of integration C[1], C[2], C[3] go away, either by 
specifying some more equations (initial conditions?) or by adding them to 
the parameter list (so that instead of solving for {k1 ... k4}, FindMinimum 
solves for {k1, ..., k4, C[1], C[2], C[3]}).  And you'll need either a fast 
machine or some serious patience.

I'm not sure, but this sounds like it should work.

-- Paul

**************************************************************************
* Paul A. Rubin                                  Phone: (517) 432-3509   *
* Department of Management                       Fax:   (517) 432-1111   *
* Eli Broad Graduate School of Management        Net:   RUBIN at MSU.EDU    *
* Michigan State University                                              *
* East Lansing, MI  48824-1122  (USA)                                    *
**************************************************************************
Mathematicians are like Frenchmen:  whenever you say something to them,
they translate it into their own language, and at once it is something
entirely different.                                    J. W. v. GOETHE


  • Prev by Date: Re: Fit in Excel Mathlink
  • Next by Date: Re: Are two expressions equal?
  • Previous by thread: Nonlinear fit to diff eqs.
  • Next by thread: Gram matrix