Matlab:求ODE系统的系数

Question

我拥有所有数据和一个由三个方程组成的 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)？

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(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')