Skip to content

Commit a501f46

Browse files
Preset for optional symbolic integration of VanderPol problem (#18)
* added the init cond with extra powwers of epsilon expansion * switched to symbolic init condition (Need to fix this) * added option for symbolic Y0 in state Validator * added preset for symbolic vanderpol integration * fixed pull request comments: conversion to double before casting to sym * fixed pull request comments: conversion to double before casting to sym * fixed pull request comments: conversion to double before casting to sym * better error msg * better error msg * Remove all symbolic dependencies * Typo fix in error message * Add comments explaining epsilon expansion Co-authored-by: Steven Roberts <sroberts994@gmail.com> Co-authored-by: Steven Roberts <steven94@vt.edu>
1 parent ca21222 commit a501f46

File tree

3 files changed

+49
-10
lines changed

3 files changed

+49
-10
lines changed

src/+otp/+utils/StructParser.m

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -39,13 +39,13 @@
3939
function isValid = validateFromString(~, x, validator)
4040
switch validator
4141
case 'numeric'
42-
isValid = isnumeric(x);
42+
isValid = isnumeric(x) || isa(x, 'sym');
4343
case 'real'
4444
isValid = isreal(x);
4545
case 'finite'
46-
isValid = all(isfinite(x));
46+
isValid = all(isfinite(x), 'all');
4747
case 'integer'
48-
isValid = isreal(x) && all(floor(x) == x);
48+
isValid = isreal(x) && all(floor(x) == x, 'all');
4949
case 'vector'
5050
isValid = isvector(x);
5151
case 'matrix'
@@ -57,13 +57,13 @@
5757
case 'scalar'
5858
isValid = isscalar(x);
5959
case 'positive'
60-
isValid = all(x > 0);
60+
isValid = all(x > 0, 'all');
6161
case 'nonnegative'
62-
isValid = all(x >= 0);
62+
isValid = all(x >= 0, 'all');
6363
case 'negative'
64-
isValid = all(x < 0);
64+
isValid = all(x < 0, 'all');
6565
case 'nonpositive'
66-
isValid = all(x <= 0);
66+
isValid = all(x <= 0, 'all');
6767
case 'row'
6868
isValid = isrow(x);
6969
case 'column'
Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
classdef ExtendedPrecision < otp.vanderpol.VanderpolProblem
2+
methods
3+
function obj = ExtendedPrecision(epsilon)
4+
if nargin < 1
5+
epsilon = 1e-6;
6+
end
7+
8+
cls = class(epsilon);
9+
10+
coeffs = [ ...
11+
% Express large fractions using numbers that fit into a double
12+
cast(199944165973, cls) / cast(3486784401, cls) * cast(7499951716, cls) / cast(4782969, cls), ...
13+
cast(-1532057748884056, cls) / cast(617673396283947, cls) * cast(4765, cls), ...
14+
cast(336712727619991, cls) / cast(22876792454961, cls) * cast(116, cls), ...
15+
cast(-231923508930412, cls) / cast(847288609443, cls), ...
16+
cast(515714746118, cls) / cast(10460353203, cls), ...
17+
cast(-3927188888, cls) / cast(387420489, cls), ...
18+
cast(34896076, cls) / cast(14348907, cls), ...
19+
cast(-1121308, cls) / cast(1594323, cls), ...
20+
cast(15266, cls) / cast(59049, cls), ...
21+
cast(-292, cls) / cast(2187, cls), ...
22+
cast(10, cls) / cast(81, cls), ...
23+
cast(-2, cls) / cast(3, cls)];
24+
z0 = 0;
25+
for c = coeffs
26+
% Use Horner's method to compute epsilon expansion
27+
z0 = epsilon * z0 + c;
28+
end
29+
30+
y0 = [2; z0];
31+
tspan = [0, 0.5];
32+
33+
params.epsilon = epsilon;
34+
obj = obj@otp.vanderpol.VanderpolProblem(tspan, y0, params);
35+
36+
end
37+
end
38+
end
39+

src/+otp/Problem.m

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -120,12 +120,12 @@ function validateNewState(obj, newTimeSpan, newY0, newParameters)
120120
% Ensures the TimeSpan, Y0, and Parameters are valid
121121
if length(newTimeSpan) ~= 2
122122
error('TimeSpan must be a vector of two times');
123-
elseif ~isnumeric(newTimeSpan)
123+
elseif ~(isnumeric(newTimeSpan) || isa(newTimeSpan, 'sym'))
124124
error('TimeSpan must be numeric');
125125
elseif ~iscolumn(newY0)
126126
error('Y0 must be a column vector');
127-
elseif ~(isnumeric(newY0) && all(isfinite(newY0)))
128-
error('Y0 must be numeric and finite');
127+
elseif ~(isnumeric(newY0) || isa(newY0, 'sym'))
128+
error('Y0 must be numeric');
129129
elseif ~(isempty(obj.ExpectedNumVars) || length(newY0) == obj.ExpectedNumVars)
130130
error('Expected Y0 to have %d components but has %d', obj.ExpectedNumVars, length(newY0));
131131
elseif ~isstruct(newParameters)

0 commit comments

Comments
 (0)