Matlab: finding coefficients of ODE system

Question

I have all the data and an ODE system of three equations which has 9 unknown coefficients (a1, a2,..., a9).

dS/dt = a1*S+a2*D+a3*F
dD/dt = a4*S+a5*D+a6*F
dF/dt = a7*S+a8*D+a9*F

t = [1 2 3 4 5]
S = [17710 18445 20298 22369 24221]
D = [1357.33 1431.92 1448.94 1388.33 1468.95]
F = [104188 104792 112097 123492 140051]

How to find these coefficients (a1,..., a9) of an ODE using Matlab?

Answer 1

I can't spend too much time on this, but basically you need to use math to reduce the equation to something more meaningful:

your equation is of the order

dx/dt = A*x

ergo the solution is

x(t-t0) = exp(A*(t-t0)) * x(t0)

Thus

exp(A*(t-t0)) = x(t-t0) * Pseudo(x(t0))

Pseudo is the Moore-Penrose Pseudo-Inverse.

EDIT: Had a second look at my solution, and I didn't calculate the pseudo-inverse properly.

Basically, Pseudo(x(t0)) = x(t0)'*inv(x(t0)*x(t0)'), as x(t0) * Pseudo(x(t0)) equals the identity matrix

Now what you need to do is assume each time step (1 to 2, 2 to 3, 3 to 4) is an experiment (therefore t-t0=1), so the solution would be to:

1- Build your pseudo inverse:

xt = [S;D;F];
xt0 = xt(:,1:4);

xInv = xt0'*inv(xt0*xt0');

2- Get exponential result

xt1 = xt(:,2:5);
expA =  xt1 * xInv;

3- Get the logarithm of the matrix:

A = logm(expA);

And since t-t0= 1, A is our solution.

And a simple proof to check

[t, y] = ode45(@(t,x) A*x,[1 5], xt(1:3,1));
plot (t,y,1:5, xt,'x')

Answer 2

You have a linear, coupled system of ordinary differential equations,

y' = Ay    with    y = [S(t);  D(t);  F(t)]

and you're trying to solve the inverse problem,

A = unknown

Interesting!

First line of attack

For given A, it is possible to solve such systems analytically (read the wiki for example).

The general solution for 3x3 design matrices A take the form

[S(t) D(t) T(t)].' = c1*V1*exp(r1*t) + c2*V2*exp(r2*t) + c3*V3*exp(r3*t)

with V and r the eigenvectors and eigenvalues of A, respectively, and c scalars that are usually determined by the problem's initial values.

Therefore, there would seem to be two steps to solve this problem:

1. Find vectors c*V and scalars r that best-fit your data
2. reconstruct A from the eigenvalues and eigenvectors.

However, going down this road is treaturous. You'd have to solve the non-linear least-squares problem for the sum-of-exponentials equation you have (using lsqcurvefit, for example). That would give you vectors c*V and scalars r. You'd then have to unravel the constants c somehow, and reconstruct the matrix A with V and r.

So, you'd have to solve for c (3 values), V (9 values), and r (3 values) to build the 3x3 matrix A (9 values) -- that seems too complicated to me.

Simpler method

There is a simpler way; use brute-force:

function test

    % find  
    [A, fval] = fminsearch(@objFcn, 10*randn(3))

end

function objVal = objFcn(A)

    % time span to be integrated over
    tspan = [1 2 3 4 5];

    % your desired data
    S = [17710    18445    20298    22369    24221   ];
    D = [1357.33  1431.92  1448.94  1388.33  1468.95 ];
    F = [104188   104792   112097   123492   140051  ];

    y_desired = [S; D; F].';

    % solve the ODE
    y0 =  y_desired(1,:);
    [~,y_real] = ode45(@(~,y) A*y, tspan, y0);

    % objective function value: sum of squared quotients
    objVal = sum((1 - y_real(:)./y_desired(:)).^2);

end

So far so good.

However, I tried both the complicated way and the brute-force approach above, but I found it very difficult to get the squared error anywhere near satisfyingly small.

The best solution I could find, after numerous attempts:

A =
    1.216731997197118e+000    2.298119167536851e-001   -2.050312097914556e-001
   -1.357306715497143e-001   -1.395572220988427e-001    2.607184719979916e-002
    5.837808840775175e+000   -2.885686207763313e+001   -6.048741083713445e-001

fval =
    3.868360951628554e-004

Which isn't bad at all :) But I would've liked a solution that was less difficult to find...