你有一个线性的、耦合的常微分方程系统,
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)
分别是V
和r
的特征向量和特征值A
,以及c
通常由问题的初始值确定的标量。
因此,似乎有两个步骤可以解决此问题:
- 查找最适合您的数据的向量
c*V
和标量r
A
从特征值和特征向量重构。
然而,沿着这条路走下去是有害的。您必须解决您拥有的指数和方程的非线性最小二乘问题(lsqcurvefit
例如,使用 )。那会给你向量c*V
和标量r
。然后你必须以某种方式解开常数,并用andc
重建矩阵。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
这一点都不错:)但我希望找到一个不太难找到的解决方案......