matlab - 结合 solve 和 dsolve 来求解具有微分和代数方程的方程系统

Question

我正在尝试求解方程系统，其中包含代数方程和微分方程。为了象征性地做到这一点，我需要将 dsolve 和 solve 结合起来(是吗？)。

考虑以下示例:我们有三个基本方程

a == b + c; % algebraic equation
diff(b,1) == 1/C1*y(t); % differential equation 1
diff(c,1) == 1/C2*y(t); % differential equation 2

求解两个微分方程，消去 int(y,0..t)，然后求解 c=f(C1,C2,a)，得到

C1*b == C2*c   or   C1*(a-c) == C2*c
c = C1/(C1+C2) * a

我怎样才能说服 Matlab 给我那个结果？这是我尝试过的:

syms a b c y C1 C2;
Eq1 = a == b + c; % algebraic equation
dEq1 = 'Db == 1/C1*y(t)'; % differential equation 1
dEq2 = 'Dc == 1/C2*y(t)'; % differential equation 2
[sol_dEq1, sol_dEq2]=dsolve(dEq1,dEq2,'b(0)==0','c(0)==0'); % this works, but no inclusion of algebraic equation
%[sol_dEq1, sol_dEq2]=dsolve(dEq1,dEq2,Eq1,'c'); % does not work
%solve(Eq1,dEq1,dEq2,'c') % does not work
%solve(Eq1,sol_dEq_C1,sol_dEq_C2,'c') % does not work

solve 和/或 dsolve 与我尝试过的方程或其解的组合没有给我有用的结果。有什么想法吗？

Answer 1

现在我假设您希望代码相当通用，所以我让它能够处理任何给定数量的方程和任何给定数量的变量，并且我没有进行手工计算。

请注意，符号工具箱的工作方式每年都会发生巨大变化，但希望这对您有用。现在可以将方程 Eq1 添加到 dSolve 的输入列表中，但这样做有两个问题:一个是 dSolve 似乎更喜欢字符输入，第二个是 dSolve 似乎没有意识到有 3 个自变量 a、b 和 c (它只看到 2 个变量，b 和 c)。

为了解决第二个问题，我对原方程微分得到一个新的微分方程，其中有三个问题:第一个是Matlab计算a关于的导数t 为 0，所以我不得不将 a 替换为 a(t) 等 b 和 c(我称 a(t) 为 a 的长版本)。第二个问题是 Matlab 使用不一致的符号，而不是将 a 的导数表示为 Da，而是将其表示为 diff(a(t), t) 因此我不得不用前者替换后者，例如 b 和 c；这给了我 Da = Db + Dc。最后一个问题是系统现在尚未确定，所以我必须获得初始值，在这里我可以解决 a(0) 但 Matlab 似乎很乐意使用 a (0) = b(0) + c(0)。

现在回到最初的第一个问题，为了解决我必须将每个符号转换回字符的问题。

这是代码

function SolveExample
syms a b c y C1 C2 t;
Eq1 = sym('a = b + c'); 
dEq1 = 'Db = 1/C1*y(t)';
dEq2 = 'Dc = 1/C2*y(t)'; 
[dEq3, initEq3] = ...
  TurnEqIntoDEq(Eq1, [a b c], t, 0);

% In the most general case Eq1 will be an array
% and thus DEq3 will be one too
dEq3_char = SymArray2CharCell(dEq3);
initEq3_char = SymArray2CharCell(initEq3);

% Below is the same as 
% dsolve(dEq1, dEq2, 'Da = Db + Dc', ...
%   'b(0)=0','c(0)=0', 'a(0) = b(0) + c(0)', 't');
[sol_dEq1, sol_dEq2, sol_dEq3] = dsolve(...
  dEq1, dEq2, dEq3_char{:}, ...
  'b(0)=0','c(0)=0', initEq3_char{:}, 't')

end

function [D_Eq, initEq] = ...
  TurnEqIntoDEq(eq, depVars, indepVar, initialVal)
% Note that eq and depVars 
% may all be vectors or scalars
% and they need not be the same size.
% eq = equations
% depVars = dependent variables
% indepVar = independent variable
% initialVal = initial value of indepVar

depVarsLong = sym(zeros(size(depVars)));
for k = 1:numel(depVars)
  % Make the variables functions
  % eg. a becomes a(t)
  % This is so that diff(a, t) does not become 0
  depVarsLong(k) = sym([char(depVars(k)) '(' ...
    char(indepVar) ')']);
end

% Next make the equation in terms of these functions
eqLong = subs(eq, depVars, depVarsLong);

% Now find the ODE corresponding to the equation
D_EqLong = diff(eqLong, indepVar);

% Now replace all the long terms like 'diff(a(t), t)'
% with short terms like 'Da'
% otherwise dSolve will not work.
% First make the short variables 'Da'
D_depVarsShort = sym(zeros(size(depVars)));
for k = 1:numel(depVars)
  D_depVarsShort(k) = sym(['D' char(depVars(k))]);
end
% Next make the long names like 'diff(a(t), t)'
D_depVarsLong = diff(depVarsLong, indepVar);
% Finally replace
D_Eq = subs(D_EqLong, D_depVarsLong, D_depVarsShort);

% Finally determine the equation
% governing the initial values
initEq = subs(eqLong, indepVar, initialVal);
end

function cc = SymArray2CharCell(sa)
cc = cell(size(sa));
for k = 1:numel(sa)
  cc{k} = char(sa(k));
end

end

一些小注意事项，我将 == 更改为 = 因为这似乎是我们 Matlab 版本之间的差异。另外，我在 dsolve 中添加了 t 作为自变量。我还假设您了解单元格、numel、线性索引等。