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