Matlab：找到ODE系统的系数

Question

我有所有数据和一个包含 9 个未知系数（a1、a2、...、a9）的三个方程的 ODE 系统。

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]

如何使用 Matlab 求 ODE 的这些系数 (a1,...,a9)？

Answer 1

我不能在这上面花太多时间，但基本上你需要使用数学来将方程简化为更有意义的东西：

你的等式是有序的

dx/dt = A*x

解决办法是

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

这样

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

Pseudo 是 Moore-Penrose Pseudo-Inverse。

编辑：再次查看我的解决方案，但我没有正确计算伪逆。

基本上，Pseudo(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')

Answer 2

你有一个线性的、耦合的常微分方程组，

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

而你正在尝试解决逆问题，

A = unknown

有趣！

第一道攻击线

对于给定的A，可以解析地求解此类系统（例如阅读the wiki）。

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，并用V 和r 重构矩阵A。

因此，您必须求解 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

这一点也不坏 :) 但我希望找到一个不太难找到的解决方案...