Combining solve and dsolve to solve equation systems with differential and algebraic equations

Question

I am trying to solve equation systems, which contain algebraic as well as differential equations. To do this symbolically I need to combine dsolve and solve (do I?).

Consider the following example: We have three base equations

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

Solving both differential equations, eliminating int(y,0..t) and then solving for c=f(C1,C2,a) yields

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

How can I convince Matlab to give me that result? Here is what I tried:

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

No combination of solve and/or dsolve with the equations or their solutions I tried gives me a useful result. Any ideas?

Answer 1

Now I assumed that you wanted the code to be rather general, so I made it to be able to work with any given number of equations and any given number of variables, and I did no calculation by hand.

Note that the way that the symbolic Toolbox works changes drastically from year to year, but hopefully this will work for you. Now one can add the equation Eq1 to the list of inputs of dSolve but there are two problems with that: One is that dSolve seems to prefer character inputs and the second is that dSolve doesn't seem to realize that there are 3 independent variables a, b, and c (it only sees 2 variables, b and c).

To solve the second problem, I differentiated the original equation to get a new differential equation, there were three problems with that: the first is Matlab evaluated the derivative of a with respect to t as 0, so I had to replace a with a(t) and such like for b and c (I called a(t) the long version of a). The second problem was Matlab used inconsistent notation, instead of representing the derivative of a as Da, it represented it as diff(a(t), t) thus I had to replace the latter with the former and such like for b and c; this gave me Da = Db + Dc. The final problem with that is that the system is now under determined, so I had to get the initial values, here I could have solved for a(0) but Matlab seemed happy with using a(0) = b(0) + c(0).

Now back to the original first problem, to solve that I had to convert every sym back into a char.

Here is the code

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

Some minor notes, I changed the == to = as that seems to be a difference between our versions of Matlab. Also I added t as the independent variable in dsolve. I also assumed that you know about cells, numel, linear indexes, ect.