2

我有所有数据和一个由三个方程组成的 ODE 系统,它有 9 个未知系数(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]

如何使用 Matlab 找到 ODE 的这些系数(a1,...,a9)?

4

2 回答 2

1

你有一个线性的、耦合的常微分方程系统,

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

你正在尝试解决问题,

A = unknown

有趣的!

第一线攻击

对于 given A,可以分析地解决此类系统(例如阅读wiki)。

3x3 设计矩阵的一般解决方案A采用以下形式

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

分别是Vr的特征向量和特征值A,以及c通常由问题的初始值确定的标量。

因此,似乎有两个步骤可以解决此问题:

  1. 查找最适合您的数据的向量c*V和标量r
  2. A从特征值和特征向量重构。

然而,沿着这条路走下去是有害的。您必须解决您拥有的指数和方程的非线性最小二乘问题(lsqcurvefit例如,使用 )。那会给你向量c*V和标量r。然后你必须以某种方式解开常数,并用andc重建矩阵。AVr

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

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

于 2013-05-15T13:48:31.583 回答
1

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

你的方程是有序的

dx/dt = A*x

因此解决方案是

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

因此

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

伪是摩尔-彭罗斯伪逆。

编辑:再看一下我的解决方案,我没有正确计算伪逆。

基本上, 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')

在此处输入图像描述

于 2013-05-15T15:29:46.090 回答