Matlab: нахождение коэффициентов системы ОДУ

Я не могу тратить на это слишком много времени, но в основном вам нужно использовать математику, чтобы сократить уравнение до чего-то более значимого:

ваше уравнение имеет порядок

dx/dt = A*x

следовательно, решение

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

Таким образом

exp(A*(t-t0)) = x(t-t0) * псевдо(x(t0))

Псевдо-это псевдообратная Мура-Пенроуза.

РЕДАКТИРОВАТЬ: я еще раз взглянул на свое решение, и я неправильно рассчитал псевдоинверсию.

По сути, Псевдо(x(t0)) = x(t0)'*inv(x(t0)*x(t0)'), как x(t0) * Pseudo(x(t0)) равно единичной матрице

Теперь вам нужно предположить, что каждый временной шаг (от 1 до 2, от 2 до 3, от 3 до 4) является экспериментом (поэтому t-t0=1), поэтому решение будет таким:

1- Создайте свою псевдоинверсию:

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

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

2- Получить экспоненциальный результат

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

3- Получить логарифм матрицы:

A = logm(expA);

А так как t-t0= 1, то A — наше решение.

И простое доказательство для проверки

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

У вас есть линейная связанная система обыкновенных дифференциальных уравнений,

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

и вы пытаетесь решить обратную задачу,

A = unknown

Интересно!

Первая линия атаки

Для заданного A такие системы можно решать аналитически (например, прочтите вики).

Общее решение для расчетных матриц 3x3 A принимает вид

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

с V и r собственными векторами и собственными значениями A соответственно и скалярами c, которые обычно определяются начальными значениями задачи.

Таким образом, казалось бы, есть два шага для решения этой проблемы:

1. Найдите векторы c*V и скаляры r, которые лучше всего соответствуют вашим данным.
2. восстановить A из собственных значений и собственных векторов.

Тем не менее, идти по этой дороге неприятно. Вам нужно будет решить нелинейную задачу наименьших квадратов для уравнения суммы экспонент, которое у вас есть (например, используя lsqcurvefit). Это даст вам векторы c*V и скаляры r. Затем вам придется каким-то образом распутать константы c и реконструировать матрицу A с помощью V и r.

Итак, вам нужно найти c (3 значения), V (9 значений) и r (3 значения), чтобы построить матрицу 3x3 A (9 значений) — мне это кажется слишком сложным.

Более простой метод

Есть более простой способ; использовать грубую силу:

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

Все идет нормально.

Тем не менее, я попробовал как сложный способ, так и подход грубой силы, описанный выше, но мне было очень трудно получить квадрат ошибки где-либо близко к удовлетворительно малой.

Лучшее решение, которое я смог найти после многочисленных попыток:

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

Что совсем неплохо :) Но мне бы хотелось, чтобы решение было проще найти...