Skip to content

Commit 0d00704

Browse files
tamaskistamaskis
committed
removed unified solve_ivp, updated other ones
1 parent 0ceee18 commit 0d00704

File tree

3 files changed

+272
-407
lines changed

3 files changed

+272
-407
lines changed

‎toolbox/ivpsolvers/solve_ivp.m‎

Lines changed: 268 additions & 67 deletions
Original file line numberDiff line numberDiff line change
@@ -1,19 +1,14 @@
11
%==========================================================================
22
%
3-
% solve_ivp Solve vector-valued and matrix-valued initial value problems
4-
% using fixed-step IVP solvers.
3+
% solve_ivp Fixed-step IVP solvers for solving vector-valued initial value
4+
% problems.
55
%
66
% [t,y] = solve_ivp(f,[t0,tf],y0,h)
77
% [t,y] = solve_ivp(f,{t0,E},y0,h)
88
% [t,y] = solve_ivp(__,method)
99
% [t,y] = solve_ivp(__,method,wb)
1010
%
11-
% [t,M] = solve_ivp(F,[t0,tf],M0,h)
12-
% [t,M] = solve_ivp(F,{t0,E},M0,h)
13-
% [t,M] = solve_ivp(__,method)
14-
% [t,M] = solve_ivp(__,method,wb)
15-
%
16-
% See also solve_ivp_matrix, solve_ivp_vector.
11+
% See also solve_ivp_matrix.
1712
%
1813
% Copyright © 2021 Tamas Kis
1914
% Last Update: 2022-09-17
@@ -28,10 +23,6 @@
2823
%
2924
%--------------------------------------------------------------------------
3025
%
31-
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
32-
% VECTOR-VALUED INITIAL VALUE PROBLEMS
33-
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
34-
%
3526
% ------
3627
% INPUT:
3728
% ------
@@ -69,67 +60,277 @@
6960
% the solution) corresponding to the nth time in "t". This convention
7061
% is chosen to match the convention used by MATLAB's ODE suite.
7162
%
72-
%--------------------------------------------------------------------------
73-
%
74-
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
75-
% MATRIX-VALUED INITIAL VALUE PROBLEMS
76-
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
77-
%
78-
% ------
79-
% INPUT:
80-
% ------
81-
% F - (1×1 function_handle) dM/dt = F(t,M) --> multivariate,
82-
% matrix-valued function (F : ℝ×ℝᵖˣʳ → ℝᵖˣʳ) defining
83-
% matrix-valued ODE
84-
% I - defines interval over which to solve the IVP, 2 options:
85-
% --> [t0,tf] - (1×2 double) initial and final times
86-
% --> {t0,E} - (1×2 cell) initial time, t₀, and function
87-
% handle for event function, E(t,M)
88-
% (E : ℝ×ℝᵖˣʳ → B)
89-
% M0 - (p×r double) initial condition, M₀ = M(t₀)
90-
% h - (1×1 double) step size
91-
% method - (OPTIONAL) (char) integration method --> 'Euler', 'RK2',
92-
% 'RK2 Heun', 'RK2 Ralston', 'RK3', 'RK3 Heun', 'RK3 Ralston',
93-
% 'SSPRK3', 'RK4', 'RK4 Ralston', 'RK4 3/8', 'AB2', 'AB3',
94-
% 'AB4', 'AB5', 'AB6', 'AB7', 'AB8', 'ABM2', 'ABM3', 'ABM4',
95-
% 'ABM5', 'ABM6', 'ABM7', 'ABM8' (defaults to 'RK4')
96-
% wb - (OPTIONAL) (1×1 logical or char) waitbar parameters (defaults
97-
% to false)
98-
% --> input as true if you want waitbar with default message
99-
% displayed
100-
% --> input as a char array storing a message if you want a
101-
% custom message displayed on the waitbar
102-
%
103-
% -------
104-
% OUTPUT:
105-
% -------
106-
% t - ((N+1)×1 double) time vector
107-
% M - (p×r×(N+1) double) solution array
108-
%
109-
% -----
110-
% NOTE:
111-
% -----
112-
% --> The nth page of "M" stores the state matrix (i.e. the solution)
113-
% corresponding to the nth time in "t".
114-
%
11563
%==========================================================================
11664
function [t,y] = solve_ivp(f,I,y0,h,method,wb)
11765

118-
% defaults optional parameters to empty vector if not input
119-
if (nargin < 5), method = []; end
120-
if (nargin < 6), wb = []; end
121-
122-
% number of columns of the state matrix
123-
r = size(y0,2);
66+
% --------------------
67+
% Time detection mode.
68+
% --------------------
12469

125-
% determines if the IVP is matrix-valued
126-
matrix_valued = (size(y0,2) > 1);
70+
if ~iscell(I)
71+
72+
% extracts initial and final times
73+
t0 = I(1);
74+
tf = I(2);
75+
76+
% makes step size negative if t0 > tf
77+
if (t0 > tf)
78+
h = -h;
79+
end
80+
81+
% defines event function
82+
if (tf > t0)
83+
E = @(t,y) t-h <= tf;
84+
else
85+
E = @(t,y) t-h >= tf;
86+
end
87+
88+
% indicates that final time is known
89+
final_time_known = true;
90+
91+
% number of subintervals between sample times
92+
N = ceil((tf-t0)/h);
93+
94+
% turns waitbar on with custom message
95+
if (nargin == 6) && ~isempty(wb) && ischar(wb)
96+
display_waitbar = true;
97+
[wb,prop] = initialize_waitbar(wb);
98+
99+
% turns waitbar on with default message
100+
elseif (nargin == 6) && ~isempty(wb) && wb
101+
display_waitbar = true;
102+
[wb,prop] = initialize_waitbar('Solving IVP...');
103+
104+
% turns waitbar off otherwise
105+
else
106+
display_waitbar = false;
107+
108+
end
109+
110+
% ---------------------
111+
% Event detection mode.
112+
% ---------------------
127113

128-
% solves IVP
129-
if matrix_valued
130-
[t,y] = solve_ivp_matrix(f,I,y0,h,method,wb);
131114
else
132-
[t,y] = solve_ivp_vector(f,I,y0,h,method,wb);
115+
116+
% extracts initial time and event function
117+
t0 = I{1};
118+
E = I{2};
119+
120+
% indicates that final time is unknown
121+
final_time_known = false;
122+
123+
end
124+
125+
% -------------------
126+
% Integration method.
127+
% -------------------
128+
129+
% defaults integration method to RK4
130+
if (nargin < 5) || isempty(method)
131+
method = 'RK4';
132+
end
133+
134+
% sets propagation function for single-step methods
135+
if strcmpi(method,'Euler')
136+
g = @(t,y) RK1_euler(f,t,y,h);
137+
elseif strcmpi(method,'RK2')
138+
g = @(t,y) RK2(f,t,y,h);
139+
elseif strcmpi(method,'RK2 Heun')
140+
g = @(t,y) RK2_heun(f,t,y,h);
141+
elseif strcmpi(method,'RK2 Ralston')
142+
g = @(t,y) RK2_ralston(f,t,y,h);
143+
elseif strcmpi(method,'RK3')
144+
g = @(t,y) RK3(f,t,y,h);
145+
elseif strcmpi(method,'RK3 Heun')
146+
g = @(t,y) RK3_heun(f,t,y,h);
147+
elseif strcmpi(method,'RK3 Ralston')
148+
g = @(t,y) RK3_ralston(f,t,y,h);
149+
elseif strcmpi(method,'SSPRK3')
150+
g = @(t,y) SSPRK3(f,t,y,h);
151+
elseif strcmpi(method,'RK4')
152+
g = @(t,y) RK4(f,t,y,h);
153+
elseif strcmpi(method,'RK4 3/8')
154+
g = @(t,y) RK4_38(f,t,y,h);
155+
elseif strcmpi(method,'RK4 Ralston')
156+
g = @(t,y) RK4_ralston(f,t,y,h);
157+
end
158+
159+
% sets propagation function for multi-step predictor methods
160+
if strcmpi(method,'AB2')
161+
g = @(t,F) AB2(f,t,F,h);
162+
elseif strcmpi(method,'AB3')
163+
g = @(t,F) AB3(f,t,F,h);
164+
elseif strcmpi(method,'AB4')
165+
g = @(t,F) AB4(f,t,F,h);
166+
elseif strcmpi(method,'AB5')
167+
g = @(t,F) AB5(f,t,F,h);
168+
elseif strcmpi(method,'AB6')
169+
g = @(t,F) AB6(f,t,F,h);
170+
elseif strcmpi(method,'AB7')
171+
g = @(t,F) AB7(f,t,F,h);
172+
elseif strcmpi(method,'AB8')
173+
g = @(t,F) AB8(f,t,F,h);
174+
end
175+
176+
% sets propagation function for multi-step predictor-corrector methods
177+
if strcmpi(method,'ABM2')
178+
g = @(t,F) ABM2(f,t,F,h);
179+
elseif strcmpi(method,'ABM3')
180+
g = @(t,F) ABM3(f,t,F,h);
181+
elseif strcmpi(method,'ABM4')
182+
g = @(t,F) ABM4(f,t,F,h);
183+
elseif strcmpi(method,'ABM5')
184+
g = @(t,F) ABM5(f,t,F,h);
185+
elseif strcmpi(method,'ABM6')
186+
g = @(t,F) ABM6(f,t,F,h);
187+
elseif strcmpi(method,'ABM7')
188+
g = @(t,F) ABM7(f,t,F,h);
189+
elseif strcmpi(method,'ABM8')
190+
g = @(t,F) ABM8(f,t,F,h);
191+
end
192+
193+
% determines if the method is a single-step method
194+
if ~strcmpi(method(1:2),'AB')
195+
single_step = true;
196+
else
197+
single_step = false;
198+
end
199+
200+
% determines order for multistep methods
201+
if ~single_step
202+
m = str2double(method(end));
203+
end
204+
205+
% --------------------------------------------------
206+
% Preallocates arrays and stores initial conditions.
207+
% --------------------------------------------------
208+
209+
% state dimension
210+
p = length(y0);
211+
212+
% preallocates time vector and solution matrix
213+
t = zeros(10000,1);
214+
y = zeros(p,10000);
215+
216+
% stores initial conditions
217+
t(1) = t0;
218+
y(:,1) = y0;
219+
220+
% ---------------------------------------
221+
% Solves IVP (using single-step methods).
222+
% ---------------------------------------
223+
224+
if single_step
225+
226+
% state vector propagation while event function is satisfied
227+
n = 1;
228+
while E(t(n),y(:,n))
229+
230+
% expands t and y if needed
231+
if (n+1) > length(t)
232+
[t,y] = expand_ivp_arrays(t,y);
233+
end
234+
235+
% propagates state vector to next sample time
236+
y(:,n+1) = g(t(n),y(:,n));
237+
238+
% increments time and loop index
239+
t(n+1) = t(n)+h;
240+
n = n+1;
241+
242+
% updates waitbar
243+
if final_time_known && display_waitbar
244+
prop = update_waitbar(n,N,wb,prop);
245+
end
246+
247+
end
248+
249+
% --------------------------------------
250+
% Solves IVP (using multi-step methods).
251+
% --------------------------------------
252+
253+
else
254+
255+
% time vector for first m+1 sample times
256+
t(1:(m+1)) = (t0:h:(t0+m*h)).';
257+
258+
% propagating state vector using RK4 for first "m" sample times
259+
for n = 1:m
260+
y(:,n+1) = RK4(f,t(n),y(:,n),h);
261+
end
262+
263+
% initializes F matrix (stores function evaluations for first m
264+
% sample times in first m columns, state vector at (m+1)th sample
265+
% time in (m+1)th column)
266+
F = zeros(p,m+1);
267+
for n = 1:m
268+
F(:,n) = f(t(n),y(:,n));
269+
end
270+
F(:,n+1) = y(:,m+1);
271+
272+
% state vector propagation while event function is satisfied
273+
n = m+1;
274+
while E(t(n),y(:,n))
275+
276+
% expands t and y if needed
277+
if (n+1) > length(t)
278+
[t,y] = expand_ivp_arrays(t,y);
279+
end
280+
281+
% updates F matrix (propagates to next sample time)
282+
F = g(t(n),F);
283+
284+
% extracts/stores state vector at next sample time
285+
y(:,n+1) = F(:,m+1);
286+
287+
% increments time and loop index
288+
t(n+1) = t(n)+h;
289+
n = n+1;
290+
291+
% updates waitbar
292+
if final_time_known && display_waitbar
293+
prop = update_waitbar(n,N,wb,prop);
294+
end
295+
296+
end
297+
298+
end
299+
300+
% -----------------
301+
% Final formatting.
302+
% -----------------
303+
304+
% trims arrays
305+
y = y(:,1:(n-1));
306+
t = t(1:(n-1));
307+
308+
% number of subintervals
309+
N = length(t)-1;
310+
311+
% linearly interpolates to find solution at desired final time
312+
if final_time_known
313+
314+
% linearly interpolates for solution at tf
315+
y(:,N+1) = y(:,N)+((y(:,N+1)-y(:,N))/(t(N+1)-t(N)))*(tf-t(N));
316+
317+
% replaces last element of "t" with tf
318+
t(N+1) = tf;
319+
320+
end
321+
322+
% deletes penultimate solution if at same time as last solution
323+
if (abs(t(N+1)-t(N)) < 1e-10)
324+
t(N) = [];
325+
y(:,N) = [];
326+
end
327+
328+
% transposes solution matrix so it is returned in "standard form"
329+
y = y.';
330+
331+
% closes waitbar
332+
if final_time_known && display_waitbar
333+
close(wb);
133334
end
134335

135336
end

‎toolbox/ivpsolvers/solve_ivp_matrix.m‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,14 @@
11
%==========================================================================
22
%
3-
% solve_ivp_matrix Solve matrix-valued initial value problems using
4-
% fixed-step IVP solvers.
3+
% solve_ivp_matrix Fixed-step IVP solvers for solving matrix-valued
4+
% initial value problems.
55
%
66
% [t,M] = solve_ivp_matrix(F,[t0,tf],M0,h)
77
% [t,M] = solve_ivp_matrix(F,{t0,E},M0,h)
88
% [t,M] = solve_ivp_matrix(__,method)
99
% [t,M] = solve_ivp_matrix(__,method,wb)
1010
%
11-
% See also solve_ivp, solve_ivp_vector.
11+
% See also solve_ivp.
1212
%
1313
% Copyright © 2021 Tamas Kis
1414
% Last Update: 2022-09-17
@@ -82,7 +82,7 @@
8282
if iscell(I), I(2) = mat2vec_E(I(2),p); end
8383

8484
% solves corresponding vector-valued IVP
85-
[t,y] = solve_ivp_vector(f,I,y0,h,method,wb);
85+
[t,y] = solve_ivp(f,I,y0,h,method,wb);
8686

8787
% transforms solution matrix for vector-valued IVP into solution array
8888
% for matrix-valued IVP

0 commit comments

Comments
 (0)