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