Skip to content

Commit d03d094

Browse files
AndreyAPopovAndreyAPopov
authored
Merge pull request #33 from ComputationalScienceLaboratory/RHS-operators
Operator overloading on RHS
2 parents 33b7f02 + 295c19a commit d03d094

File tree

3 files changed

+172
-2
lines changed

3 files changed

+172
-2
lines changed

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

+3-1
Original file line numberDiff line numberDiff line change
@@ -3,7 +3,9 @@
33
function obj = Canonical()
44
tspan = [0.5, 2];
55
y0 = [sinh(0.5); tanh(0.5)];
6-
obj = [email protected](tspan, y0, struct());
6+
params = otp.trigonometricdae.TrigonometricDAEParameters;
7+
obj = [email protected](tspan, ...
8+
y0, params);
79
end
810
end
911
end

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

+4
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,4 @@
1+
classdef TrigonometricDAEParameters
2+
%TRIGONOMETRICDAEPARAMETERS
3+
end
4+

‎src/+otp/RHS.m

+165-1
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
classdef RHS
1+
classdef RHS
22
properties (SetAccess = private)
33
F
44

@@ -22,6 +22,13 @@
2222
OnEvent
2323
end
2424

25+
properties (Dependent)
26+
JacobianMatrix
27+
JacobianFunction
28+
MassMatrix
29+
MassFunction
30+
end
31+
2532
methods
2633
function obj = RHS(F, varargin)
2734
obj.F = F;
@@ -34,6 +41,69 @@
3441
end
3542
end
3643

44+
function mat = get.JacobianMatrix(obj)
45+
mat = obj.prop2Matrix(obj.Jacobian);
46+
end
47+
48+
function fun = get.JacobianFunction(obj)
49+
fun = obj.prop2Function(obj.Jacobian);
50+
end
51+
52+
function mat = get.MassMatrix(obj)
53+
mat = obj.prop2Matrix(obj.Mass);
54+
end
55+
56+
function fun = get.MassFunction(obj)
57+
fun = obj.prop2Function(obj.Mass);
58+
end
59+
60+
function obj = uplus(obj)
61+
end
62+
63+
function newRHS = uminus(obj)
64+
newRHS = mtimes(-1, obj);
65+
end
66+
67+
function newRHS = plus(obj1, obj2)
68+
newRHS = applyOp(obj1, obj2, @plus, @(r, ~) r, @(~, r) r, @plus);
69+
end
70+
71+
function newRHS = minus(obj1, obj2)
72+
newRHS = applyOp(obj1, obj2, @minus, @(r, ~) r, @(~, r) -r, @minus);
73+
end
74+
75+
function newRHS = mtimes(obj1, obj2)
76+
newRHS = applyOp(obj1, obj2, @mtimes, @mtimes, @mtimes, []);
77+
end
78+
79+
function newRHS = times(obj1, obj2)
80+
newRHS = applyOp(obj1, obj2, @times, @times, @times, []);
81+
end
82+
83+
function newRHS = rdivide(obj1, obj2)
84+
newRHS = applyOp(obj1, obj2, @rdivide, @rdivide, @rdivide, []);
85+
end
86+
87+
function newRHS = ldivide(obj1, obj2)
88+
newRHS = applyOp(obj1, obj2, @ldivide, @ldivide, @ldivide, []);
89+
end
90+
91+
function newRHS = mrdivide(obj1, obj2)
92+
newRHS = applyOp(obj1, obj2, @mrdivide, @mrdivide, @mrdivide, []);
93+
end
94+
95+
function newRHS = mldivide(obj1, obj2)
96+
newRHS = applyOp(obj1, obj2, @mldivide, @mldivide, @mldivide, []);
97+
end
98+
99+
function newRHS = power(obj1, obj2)
100+
newRHS = applyOp(obj1, obj2, @power, [], [], []);
101+
end
102+
103+
function newRHS = mpower(obj1, obj2)
104+
newRHS = applyOp(obj1, obj2, @mpower, [], [], []);
105+
end
106+
37107
function opts = odeset(obj, varargin)
38108
opts = odeset( ...
39109
'Events', obj.Events, ...
@@ -48,5 +118,99 @@
48118
'Vectorized', obj.Vectorized, ...
49119
varargin{:});
50120
end
121+
122+
end
123+
124+
methods (Access = private)
125+
function mat = prop2Matrix(~, p)
126+
if isa(p, 'function_handle')
127+
mat = [];
128+
else
129+
mat = p;
130+
end
131+
end
132+
133+
function fun = prop2Function(~, p)
134+
if isa(p, 'function_handle') || isempty(p)
135+
fun = p;
136+
else
137+
fun = @(varargin) p;
138+
end
139+
end
140+
141+
function newRHS = applyOp(obj1, obj2, op, dOpLeft, dOpRight, dOpBoth)
142+
if isa(obj1, 'function_handle')
143+
obj1 = otp.RHS(obj1);
144+
elseif isa(obj2, 'function_handle')
145+
obj2 = otp.RHS(obj2);
146+
end
147+
148+
if isa(obj1, 'otp.RHS')
149+
primaryRHS = obj1;
150+
if isa(obj2, 'otp.RHS')
151+
f = otp.RHS.mergeProp(obj1.F, obj2.F, op);
152+
merge = @(p) otp.RHS.mergeProp(obj1.(p), obj2.(p), dOpBoth);
153+
154+
if strcmp(obj1.Vectorized, obj2.Vectorized)
155+
vectorized = obj1.Vectorized;
156+
end
157+
else
158+
f = otp.RHS.mergeProp(obj1.F, obj2, op);
159+
merge = @(p) otp.RHS.mergeProp(obj1.(p), obj2, dOpLeft);
160+
vectorized = obj1.Vectorized;
161+
end
162+
else
163+
primaryRHS = obj2;
164+
f = otp.RHS.mergeProp(obj1, obj2.F, op);
165+
merge = @(p) otp.RHS.mergeProp(obj1, obj2.(p), dOpRight);
166+
vectorized = obj2.Vectorized;
167+
end
168+
169+
% Events and NonNegative practically cannot be supported and are
170+
% always unset.
171+
172+
% JPattern is problematic to compute for division operators due to
173+
% singular patterns
174+
175+
% Mass matrices introduce several difficulties. When singular, it
176+
% makes it infeasible to update InitialSlope, and therefore, it is
177+
% always unset. To avoid issues with two RHS' having different mass
178+
% matrices, only the primary RHS is used.
179+
180+
newRHS = otp.RHS(f, ...
181+
'Mass', primaryRHS.Mass, ...
182+
'MassSingular', primaryRHS.MassSingular, ...
183+
'MStateDependence', primaryRHS.MStateDependence, ...
184+
'MvPattern', primaryRHS.MvPattern, ...
185+
'Jacobian', merge('Jacobian'), ...
186+
'JacobianVectorProduct', merge('JacobianVectorProduct'), ...
187+
'JacobianAdjointVectorProduct', ...
188+
merge('JacobianAdjointVectorProduct'), ...
189+
'PartialDerivativeParameters', ...
190+
merge('PartialDerivativeParameters'), ...
191+
'PartialDerivativeTime', merge('PartialDerivativeTime'), ...
192+
'HessianVectorProduct', merge('HessianVectorProduct'), ...
193+
'HessianAdjointVectorProduct', ...
194+
merge('HessianAdjointVectorProduct'), ...
195+
'Vectorized', vectorized);
196+
end
197+
end
198+
199+
methods (Static, Access = private)
200+
function p = mergeProp(p1, p2, op)
201+
if isempty(p1) || isempty(p2) || isempty(op)
202+
p = [];
203+
elseif isa(p1, 'function_handle')
204+
if isa(p2, 'function_handle')
205+
p = @(varargin) op(p1(varargin{:}), p2(varargin{:}));
206+
else
207+
p = @(varargin) op(p1(varargin{:}), p2);
208+
end
209+
elseif isa(p2, 'function_handle')
210+
p = @(varargin) op(p1, p2(varargin{:}));
211+
else
212+
p = op(p1, p2);
213+
end
214+
end
51215
end
52216
end

0 commit comments

Comments
 (0)