Info

此问题已关闭。 请重新打开它进行编辑或回答。

Did matlab really fail to solve this system of 4 DAEs or did I make a mistake somewhere?

1 次查看（过去 30 天）
Dursman Mchabe
Dursman Mchabe 2018-9-3
关闭： MATLAB Answer Bot 2021-8-20
Hi everyone, I'm trying to analytically solve a system of 4 DAEs, however all outputs are blank. Please see the code below:
function WaterLiquidFilmProfilesAnalytical
syms a b c e f t A B C E F G H K
dEq1 = 'A * D2a + B * D2b + C * D2c = 0';
dEq2 = 'E * D2e - B * D2b - 2 * C * D2c - F * D2f = 0';
Eq1 = ('G * a = e * b');
Eq2 = str2sym(Eq1);
[dEq3, initEq3] = TurnEqIntoDEq(Eq2, [G a e b], t, 0);
dEq3_char = SymArray2CharCell(dEq3);
initEq3_char = SymArray2CharCell(initEq3);
Eq4 = ('H * b = e * c');
Eq5 = str2sym(Eq4);
[dEq4, initEq4] = TurnEqIntoDEq(Eq5, [H b e c], t, 0);
dEq4_char = SymArray2CharCell(dEq4);
initEq4_char = SymArray2CharCell(initEq4);
Eq6 = ('K = e * f');
Eq7 = str2sym(Eq6);
[dEq5, initEq5] = TurnEqIntoDEq(Eq7, [K e f], t, 0);
dEq5_char = SymArray2CharCell(dEq5);
initEq5_char = SymArray2CharCell(initEq5);
%[sol_dEq1, sol_dEq2, sol_dEq3, sol_dEq4, sol_dEq5] = dsolve(dEq1, dEq2, dEq3_char{:},'a(0)=0','e(0)=0','b(0)=0', initEq3_char{:}, 't', dEq4_char{:},'b(0)=0','e(0)=0','c(0)=0', initEq4_char{:}, 't', dEq5_char{:},'d(0)=0','f(0)=0', initEq5_char{:}, 't')
%[sol_dEq1, sol_dEq2, sol_dEq3, sol_dEq4, sol_dEq5] = dsolve(dEq1, dEq2, dEq3_char{:}, dEq4_char{:}, dEq5_char{:},'a(0)=0','b(0)=0','c(0)=0','e(0)=0','f(0)=0',initEq3_char{:},initEq4_char{:}, initEq5_char{:},'t')
[a, b, c, d, e] = dsolve(dEq1, dEq2, dEq3_char{:}, dEq4_char{:}, dEq5_char{:},'a(0)=0','b(0)=0','c(0)=0','e(0)=0','f(0)=0',initEq3_char{:},initEq4_char{:}, initEq5_char{:},'t')
end
function [D_Eq, initEq] = TurnEqIntoDEq(eq, depVars, indepVar, initialVal)
% eq = equations
% depVars = dependent variables
% indepVar = independent variable
% initialVal = initial value of indepVar
depVarsLong = sym(zeros(size(depVars)));
for k = 1:numel(depVars)
% Making the variables functions
% eg. a becomes a(t)
% This is so that diff(a, t) does not become 0
depVarsLong(k) = str2sym([char(depVars(k)) '(' char(indepVar) ')']);
end
% Next making 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'
D_depVarsShort = sym(zeros(size(depVars)));
for k = 1:numel(depVars)
D_depVarsShort(k) = str2sym(['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
Did matlab really fail or did I make a mistake?
I also tried to solve these DAEs numerically,using:
syms a(x) b(x) c(x) d(x) e(x) A B C D E FF G H
eqn1 = A * diff(a(x),x,2) + B * diff(b(x),x,2) + C * diff(c(x),x,2) == 0;
eqn2 = D * diff(d(x),x,2) - B * diff(b(x),x,2) - 2 * C * diff(c(x),x,2) - E * diff(e(x),x,2) == 0;
eqn3 = FF == (d(x) * b(x)) / a(x);
eqn4 = G == (d(x) * c(x)) / b(x);
eqn5 = H == d(x) * e(x);
eqns = [eqn1 eqn2 eqn3 eqn4 eqn5];
vars = [a(x); b(x); c(x); d(x); e(x)];
origVars = length(vars);
M = incidenceMatrix(eqns, vars);
[eqns, vars] = reduceDifferentialOrder(eqns, vars);
isLowIndexDAE(eqns,vars);
f = daeFunction(eqns,vars, A, B, C, D, E, FF, G, H);
A = 2.89e-9;
B = 2.35e-9;
C = 1.69e-9;
D = 1.6e-8;
E = 9.25e-9;
FF = 6.24;
G = 5.68e-5;
H = 5.3e-8;
F = @(t, Y, YP) f(t, Y, YP, A, B, C, D, E, FF, G, H);
vars;
y0est = [1.2650e4; 2.2748e4; 0.3724; 3.47; 1.5274e-8];
yp0est = zeros(5,1);
opt = odeset('RelTol', 10^(-3), 'AbsTol' , 10^(-3));
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);
[tSol,ySol] = ode15i(F, [0, 7.21e-05], y0, yp0, opt);
for k = 1:origVars
S{k} = char(vars(k));
end
plot(xSol,ySol(:,1),'r:',xSol,ySol(:,2),'k-.',xSol,ySol(:,3),'b--', xSol,ySol(:,4),'c-',xSol,ySol(:,5),'m-')
ylabel('Conc (mol/m^{3}')
xlabel ('Time (sec)')
legend('SO_2','HSO_{3}^{-}', 'SO_{3}^{2-}', 'H^+','OH^-','location','Best')
However, I get the error message:
Index in position 1 exceeds array bounds (must not exceed 5).
Error in
symengine>@(x,in2,in3,param1,param2,param3,param4,param5,param6,param7,param8)[in3(6,:).*param1+in3(7,:).*param2+in3(8,:).*param3;-in3(7,:).*param2-in3(8,:).*param3.*2.0+in3(9,:).*param4-in3(10,:).*param5;param6-(in2(2,:).*in2(4,:))./in2(1,:);param7-(in2(3,:).*in2(4,:))./in2(2,:);param8-in2(4,:).*in2(5,:);in2(6,:)-in3(1,:);in2(7,:)-in3(2,:);in2(8,:)-in3(3,:);in2(9,:)-in3(4,:);in2(10,:)-in3(5,:)]
Error in WaterProfiles>@(t,Y,YP)f(t,Y,YP,A,B,C,D,E,FF,G,H)
Error in decic (line 66)
res = feval(odefun,t0,y0,yp0,varargin{:});
Error in WaterProfiles (line 45)
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);
Debugging points me to line 45:
[y0, yp0] = decic(F, 0, y0est, [], yp0est, [], opt);
And in my understanding, the index of y0est does not exceed 5:
y0est = [1.2650e4; 2.2748e4; 0.3724; 3.47; 1.5274e-8];
  7 个评论
显示 5更早的评论
Dursman Mchabe
Dursman Mchabe 2018-9-4
My wish was to get an analytical solution. For consistency with other sections of my work. Because, if I use a numerical solution, I will be forced to solve a three-point bvp problem. Which is very very difficult.

此问题已关闭。

回答（0 个）

此问题已关闭。

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by