Skip to content

Commit 81ef2ab

Browse files
elswitArash SarsharSteven-Roberts
authored
2 variable index-1 DAE problem proposed by Ascher et al. (#5)
* 2 variable index-1 DAE problem proposed by Ascher et al. * Code cleanup based on PR #5 and add more mass matrix properties to Rhs * Add beta parameter to preset * Update mass matrix properties for Rhs Co-authored-by: Arash Sarshar <arash@Arashs-MacBook-Pro.local> Co-authored-by: Steven Roberts <sroberts994@gmail.com>
1 parent 2fc2f11 commit 81ef2ab

File tree

10 files changed

+98
-31
lines changed

10 files changed

+98
-31
lines changed

‎src/+otp/+ascherlineardae/+presets/Canonical.m‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
classdef Canonical < otp.ascherlineardae.AscherLinearDAEProblem
2+
methods
3+
function obj = Canonical(beta)
4+
tspan = [0.0; 1];
5+
6+
if nargin < 1
7+
beta = 0.5;
8+
end
9+
10+
params.beta = beta;
11+
12+
y0 = [1; beta];
13+
14+
obj = obj@otp.ascherlineardae.AscherLinearDAEProblem(tspan, y0, params);
15+
end
16+
end
17+
end

‎src/+otp/+ascherlineardae/AscherLinearDAEProblem.m‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,36 @@
1+
classdef AscherLinearDAEProblem < otp.Problem
2+
%Index-1 DAE problem
3+
% Problem description and fortran implementation can be found in
4+
% Ascher, Uri. "On symmetric schemes and differential-algebraic equations."
5+
% SIAM journal on scientific and statistical computing 10.5 (1989): 937-949.
6+
7+
methods
8+
function obj = AscherLinearDAEProblem(timeSpan, y0, parameters)
9+
obj@otp.Problem('Ascher Linear DAE', 2, timeSpan, y0, parameters);
10+
end
11+
end
12+
13+
methods (Access = protected)
14+
function onSettingsChanged(obj)
15+
beta = obj.Parameters.beta;
16+
17+
obj.Rhs = otp.Rhs(@(t, y) otp.ascherlineardae.f(t, y, beta), ...
18+
otp.Rhs.FieldNames.Jacobian, @(t, y) otp.ascherlineardae.jac(t, y, beta), ...
19+
otp.Rhs.FieldNames.Mass, @(t) otp.ascherlineardae.mass(t, [], beta), ...
20+
otp.Rhs.FieldNames.MStateDependence, 'none', ...
21+
otp.Rhs.FieldNames.MassSingular, 'yes');
22+
end
23+
24+
function validateNewState(obj, newTimeSpan, newY0, newParameters)
25+
validateNewState@otp.Problem(obj, newTimeSpan, newY0, newParameters)
26+
27+
otp.utils.StructParser(newParameters) ...
28+
.checkField('beta', 'scalar', 'real', 'finite');
29+
end
30+
31+
function sol = internalSolve(obj, varargin)
32+
sol = internalSolve@otp.Problem(obj, 'Method', @ode15s, varargin{:});
33+
end
34+
end
35+
end
36+

‎src/+otp/+ascherlineardae/f.m‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
function dy = f(t, y, beta)
2+
3+
dy = [-1, 1 + t; beta, -1 - beta * t] * y + [0; sin(t)];
4+
5+
end

‎src/+otp/+ascherlineardae/jac.m‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
function J = jac(t, ~, beta)
2+
3+
J = [-1, 1 + t; beta, -1 - beta * t];
4+
5+
end

‎src/+otp/+ascherlineardae/mass.m‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
function M = mass(t, ~, ~)
2+
3+
M = [1, -t; 0, 0];
4+
5+
end

‎src/+otp/+transistoramplifier/TransistorAmplifierProblem.m‎

Lines changed: 8 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -16,23 +16,24 @@ function onSettingsChanged(obj)
1616

1717
obj.Rhs = otp.Rhs(@(t, y) otp.transistoramplifier.f(t, y, C, R, Ub, UF, alpha, beta), ...
1818
otp.Rhs.FieldNames.Jacobian, @(t, y) otp.transistoramplifier.jac(t, y, C, R, Ub, UF, alpha, beta), ...
19-
otp.Rhs.FieldNames.MassMatrix, otp.transistoramplifier.mass([], [], C, R, Ub, UF, alpha, beta));
19+
otp.Rhs.FieldNames.Mass, otp.transistoramplifier.mass([], [], C, R, Ub, UF, alpha, beta), ...
20+
otp.Rhs.FieldNames.MassSingular, 'yes');
2021
end
2122

2223
function validateNewState(obj, newTimeSpan, newY0, newParameters)
2324
validateNewState@otp.Problem(obj, newTimeSpan, newY0, newParameters)
2425

2526
otp.utils.StructParser(newParameters) ...
26-
.checkField('C', 'real', 'finite', 'nonnegative') ...
27+
.checkField('C', 'real', 'finite', 'nonnegative') ...
2728
.checkField('R', 'real', 'finite', 'nonnegative') ...
28-
.checkField('Ub', 'scalar', 'real', 'finite', 'nonnegative') ...
29-
.checkField('UF', 'scalar', 'real', 'finite', 'nonnegative') ...
30-
.checkField('alpha', 'scalar', 'real', 'finite', 'nonnegative') ...
31-
.checkField('beta', 'scalar', 'real', 'finite', 'nonnegative');
29+
.checkField('Ub', 'scalar', 'real', 'finite', 'nonnegative') ...
30+
.checkField('UF', 'scalar', 'real', 'finite', 'nonnegative') ...
31+
.checkField('alpha', 'scalar', 'real', 'finite', 'nonnegative') ...
32+
.checkField('beta', 'scalar', 'real', 'finite', 'nonnegative');
3233
end
3334

3435
function sol = internalSolve(obj, varargin)
35-
sol = internalSolve@otp.Problem(obj, 'Method', @ode15s, 'AbsTol', 1e-10, varargin{:});
36+
sol = internalSolve@otp.Problem(obj, 'Method', @ode15s, varargin{:});
3637
end
3738
end
3839
end

‎src/+otp/+trigonometricdae/TrigonometricDAEProblem.m‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -6,10 +6,11 @@
66
end
77

88
methods (Access = protected)
9-
function onSettingsChanged(obj)
9+
function onSettingsChanged(obj)
1010
obj.Rhs = otp.Rhs(@(t, y) otp.trigonometricdae.f(t, y), ...
1111
otp.Rhs.FieldNames.Jacobian, @(t, y) otp.trigonometricdae.jac(t, y), ...
12-
otp.Rhs.FieldNames.MassMatrix, otp.trigonometricdae.mass([], []));
12+
otp.Rhs.FieldNames.Mass, otp.trigonometricdae.mass([], []), ...
13+
otp.Rhs.FieldNames.MassSingular, 'yes');
1314
end
1415

1516
function sol = internalSolve(obj, varargin)

‎src/+otp/+zlakinetics/ZLAKineticsProblem.m‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,8 @@ function onSettingsChanged(obj)
2020

2121
obj.Rhs = otp.Rhs(@(t, y) otp.zlakinetics.f(t, y, k, K, klA, Ks, pCO2, H), ...
2222
otp.Rhs.FieldNames.Jacobian, @(t, y) otp.zlakinetics.jac(t, y, k, K, klA, Ks, pCO2, H), ...
23-
otp.Rhs.FieldNames.MassMatrix, otp.zlakinetics.mass([], [], k, K, klA, Ks, pCO2, H));
23+
otp.Rhs.FieldNames.Mass, otp.zlakinetics.mass([], [], k, K, klA, Ks, pCO2, H), ...
24+
otp.Rhs.FieldNames.MassSingular, 'yes');
2425
end
2526

2627
function validateNewState(obj, newTimeSpan, newY0, newParameters)

‎src/+otp/Problem.m‎

Lines changed: 2 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -257,7 +257,7 @@ function validateNewState(obj, newTimeSpan, newY0, newParameters)
257257
p.addParameter('Method', @ode45);
258258
p.parse(varargin{:});
259259

260-
options = otp.Problem.odeset(obj, p.Unmatched);
260+
options = obj.Rhs.odeset(p.Unmatched);
261261

262262
sol = p.Results.Method(obj.Rhs.F, obj.TimeSpan, obj.Y0, options);
263263

@@ -273,7 +273,7 @@ function validateNewState(obj, newTimeSpan, newY0, newParameters)
273273
return;
274274
end
275275

276-
options = otp.Problem.odeset(problem, options);
276+
options = problem.Rhs.odeset(options);
277277
sol = odextend(sol, problem.Rhs.F, problem.TimeSpan(end), problem.Y0, options);
278278
end
279279
end
@@ -310,24 +310,6 @@ function validateNewState(obj, newTimeSpan, newY0, newParameters)
310310
end
311311
end
312312

313-
methods (Static, Access = private)
314-
function newOptions = odeset(problem, options)
315-
newOptions = odeset(options);
316-
317-
if isprop(problem.Rhs, 'Jacobian')
318-
newOptions.Jacobian = problem.Rhs.Jacobian;
319-
end
320-
321-
if isprop(problem.Rhs, 'MassMatrix')
322-
newOptions.Mass = problem.Rhs.MassMatrix;
323-
end
324-
325-
if isprop(problem.Rhs, 'Events')
326-
newOptions.Events = problem.Rhs.Events;
327-
end
328-
end
329-
end
330-
331313
methods (Abstract, Access = protected)
332314
% This method is called when either TimeSpan, Y0, or parameters are changed. It should update F and other properties such as a Jacobian to reflect the changes.
333315
onSettingsChanged(obj);

‎src/+otp/Rhs.m‎

Lines changed: 15 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,9 @@
1212
'JacobianTime', 'JacobianTime', ...
1313
'HessianVectorProduct', 'HessianVectorProduct', ...
1414
'HessianAdjointVectorProduct', 'HessianAdjointVectorProduct', ...
15-
'MassMatrix', 'MassMatrix',...
15+
'Mass', 'Mass',...
16+
'MStateDependence', 'MStateDependence', ...
17+
'MassSingular', 'MassSingular', ...
1618
'JPattern', 'JPattern', ...
1719
'Events', 'Events', ...
1820
'OnEvent', 'OnEvent');
@@ -41,5 +43,17 @@
4143
prop.SetAccess = 'private';
4244
end
4345
end
46+
47+
function s = struct(obj)
48+
props = properties(obj);
49+
s = struct();
50+
for i = 1:length(props)
51+
s.(props{i}) = obj.(props{i});
52+
end
53+
end
54+
55+
function opts = odeset(obj, varargin)
56+
opts = odeset(struct(obj), varargin{:});
57+
end
4458
end
4559
end

0 commit comments

Comments
 (0)