我想求解一个参数未知的ODE系统k1,k2,k3：
dC1/dx=-k1*C1
dC2/dx=k1*C1-k2*C2
dC3/dx=k2*C2-k3*C3
我有一组实验结果，C1,C2,C3 的值在x=0（入口点）和x=1（终点）。
我在x=0 和x=1 之间没有任何数据来使用像ode45() 或ode23() 这样的ODE 函数来解决它，然后使用优化函数。如何在 MATLAB 中解决这个问题？
【问题讨论】：
标签： matlab optimization model curve-fitting ode
您可以尝试使用边界值求解器bvp4c（或bvp5c）。可以满足的边界条件的数量是 ODE 的维数和未知参数的数量之和，因此您希望将参数固定为 6 个已知值的 6 个可能条件。也就是说，传递给求解器的函数是
function dCdx = odesys(x,C,k)
dCdx = zeros_like(C)
dCdx(1) = -k(1)*C(1)
dCdx(2) = k(1)*C(1)-k(2)*C(2)
dCdx(3) = k(2)*C(2)-k(3)*C(3)
function bc = boundary(C0, CT, k)
bc = [ C0(1)-c01, C0(2)-c02, C0(3)-c03, CT(1)-cT1, CT(2)-cT2, CT(3)-cT3 ]
作为初始猜测，您可以使用已知值之间的一些线性插值。
关于提供边界值的修改系统
dC1/dx=-k1*C1/(1+k1*C1) ;
dC2/dx=k1*C1/(1+k1*C1)-k2*C2 ;
dC3/dx=k2*C2-k3*C3 ;
在脚本中使用类似的 python 求解器（但将 x1 更改为 100）
c0A, c0B, c0C = 293.3 , 2.1414, 3.6884
c1A, c1B, c1C = 208.09, 33.823, 78.561
x0, x1 = 0.0, 100.0
def odesys(x,C,k):
dCdx = zeros_like(C)
dCdx[0] = -k[0]*C[0]/(1+k[0]*C[0]) ;
dCdx[1] = k[0]*C[0]/(1+k[0]*C[0])-k[1]*C[1]
dCdx[2] = k[1]*C[1]-k[2]*C[2]
return dCdx
def bc(C0, C1, k):
return [ C0[0]-c0A, C0[1]-c0B, C0[2]-c0C, C1[0]-c1A, C1[1]-c1B, C1[2]-c1C ]
x_init = linspace(x0,x1,21)
s = (x_init-x0)/(x1-x0)
C_init = [ c0A+s*(c1A-c0A), c0B+s*(c1B-c0B), c0C+s*(c1C-c0C)]
k = [1.0e-2, 1.0, 1.0]
res = solve_bvp(odesys, bc, x_init, C_init, k)
print res.message, res.p
它将参数[ k1, k2, k3]的结果打印为
The algorithm converged to the desired accuracy.
[ 0.02319266 0.02248122 -0.00678455]
第三个参数是非物理负数，暗示模型或某些参数不正确。此外，绘制解决方案：
x_sol=np.linspace(x0,x1,301)
C_sol=res.sol(x_sol)
for k in range(3): plt.subplot(1,3,k+1); plt.plot(x_sol, C_sol[k]); plt.grid()
plt.show()
似乎显示了正确的行为。使用这些参数与另一种积分方法进行前向积分可以确认解决方案的正确性
C_int = odeint( lambda C,t: odesys(t,C,res.p), res.sol(0), x_sol)
for k in range(3): plt.subplot(1,3,k+1); plt.plot(x_sol, C_int[:,k]); plt.grid()
plt.show()
【讨论】：
k 并返回端点到给定端点的距离之和。然后您可以使用通用优化求解器进行调用。
这是一个由 3 个常微分方程组成的系统问题。在x=0 和x=1 处的C1、C2 和C3 的值是初始条件。您可以使用dsolve 解析解决它。在我的示例中，我假设初始条件为C1(0)=1; C2(0)=1; C3(1)=1;。您可以根据自己的数据进行修改。
syms C1(x) C2(x) C3(x) k1 k2 k3
% Define the equations using == and represent differentiation using the diff function.
ode1 = diff(C1) == -k1*C1;
ode2 = diff(C2) == k1*C1-k2*C2;
ode3 = diff(C3) == k2*C2-k3*C3;
odes = [ode1; ode2; ode3]
% Define the initial conditions using ==
cond1 = C1(0) == 1;
cond2 = C2(0) == 1;
cond3 = C3(1) == 1;
conds = [cond1; cond2; cond3];
% dsolve function finds values for the constants that satisfy these conditions.
[C1Sol(x) C2Sol(x) C3Sol(x)] = dsolve(odes, conds)
输出
odes(x) =
D(C1)(x) == -k1*C1(x)
D(C2)(x) == k1*C1(x) - k2*C2(x)
D(C3)(x) == k2*C2(x) - k3*C3(x)
C1Sol(x) =
exp(-k1*x)
C2Sol(x) =
-(exp(-k1*x)*exp(-k2*x)*((k1^2*k2*exp(k2*x))/((k1 - k2)*(k1 - k3)) - (k2^2*exp(k1*x)*(2*k1 - k2))/((k1 - k2)*(k2 - k3)) - (k1*k2*k3*exp(k2*x))/((k1 - k2)*(k1 - k3)) + (k2*k3*exp(k1*x)*(2*k1 - k2))/((k1 - k2)*(k2 - k3))))/k2
C3Sol(x) =
-exp(-k1*x)*exp(-k2*x)*exp(-k3*x)*((k2*exp(k1*x)*exp(k3*x)*(2*k1 - k2))/((k1 - k2)*(k2 - k3)) - (k1*k2*exp(k2*x)*exp(k3*x))/((k1 - k2)*(k1 - k3)) + (exp(k1*x)*exp(k2*x)*exp(k3)*(k1*k2^2 - k1^2*k2 - k1*k3^2 + k1^2*k3 + k2*k3^2 - k2^2*k3 + k1*k2^2*exp(-k1) + k1*k2^2*exp(-k2) - 2*k1^2*k2*exp(-k2) - k2^2*k3*exp(-k2) - k1*k2*k3*exp(-k1) + 2*k1*k2*k3*exp(-k2)))/((k1 - k2)*(k1 - k3)*(k2 - k3)))
【讨论】：