Skip to content

Commit 22e6fc9

Browse files
tamaskistamaskis
committed
single unified ode solver function
1 parent e11b51c commit 22e6fc9

File tree

5 files changed

+892
-105
lines changed

5 files changed

+892
-105
lines changed

‎toolbox/odesolvers/ode.asv‎

Lines changed: 287 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,287 @@
1+
%==========================================================================
2+
%
3+
% odeRK Fixed-step ODE solver using the Runge-Kutta methods.
4+
%
5+
% [t,y] = odeRK(f,[t0,tf],y0,h)
6+
% [t,y] = odeRK(f,{t0,C},y0,h)
7+
% [t,y] = odeRK(__,method)
8+
% [t,y] = odeRK(__,method,wb)
9+
%
10+
% See also odeAB, odeABM.
11+
%
12+
% Copyright © 2021 Tamas Kis
13+
% Last Update: 2022-06-04
14+
% Website: https://tamaskis.github.io
15+
% Contact: tamas.a.kis@outlook.com
16+
%
17+
% TOOLBOX DOCUMENTATION:
18+
% https://tamaskis.github.io/ODE_Solver_Toolbox-MATLAB/
19+
%
20+
% TECHNICAL DOCUMENTATION:
21+
% https://tamaskis.github.io/documentation/Fixed_Step_ODE_Solvers.pdf
22+
%
23+
%--------------------------------------------------------------------------
24+
%
25+
% ------
26+
% INPUT:
27+
% ------
28+
% f - (1×1 function_handle) dy/dt = f(t,y) --> multivariate,
29+
% vector-valued function (f : ℝ×ℝᵖ → ℝᵖ) defining ODE
30+
% I - defines interval over which to solve the ODE, 2 options:
31+
% --> [t0,tf] - (1×2 double) initial and final times
32+
% --> {t0,C} - (1×2 cell) initial time, t₀, and function
33+
% handle for condition function, C(t,y)
34+
% (C : ℝ×ℝᵖ → 𝔹)
35+
% y0 - (p×1 double) initial condition, y₀ = y(t₀)
36+
% h - (1×1 double) step size
37+
% method - (char) (OPTIONAL) Runge-Kutta method --> 'Euler', 'RK2',
38+
% 'RK2 Heun', 'RK2 Ralston', 'RK3', 'RK3 Heun', 'RK3 Ralston',
39+
% 'SSPRK3', 'RK4', 'RK4 Ralston', 'RK4 3/8' (defaults to 'RK4')
40+
% wb - (1×1 logical or char) (OPTIONAL) waitbar parameters
41+
% --> input as "true" if you want waitbar with default
42+
% message displayed
43+
% --> input as a char array storing a message if you want a
44+
% custom message displayed on the waitbar
45+
%
46+
% -------
47+
% OUTPUT:
48+
% -------
49+
% t - ((N+1)×1 double) time vector
50+
% y - ((N+1)×p double) matrix storing time history of state vector
51+
%
52+
% -----
53+
% NOTE:
54+
% -----
55+
% --> p = dimension of state vector (for the scalar case, p = 1)
56+
% --> N+1 = length of time vector
57+
% --> The ith row of "y" is the TRANSPOSE of the state vector (i.e. the
58+
% solution corresponding to the ith time in "t"). This convention is
59+
% chosen to match the convention used by MATLAB's ODE suite.
60+
%
61+
%==========================================================================
62+
function [t,y] = ode(f,I,y0,h,method,wb)
63+
64+
% --------------------
65+
% Time detection mode.
66+
% --------------------
67+
68+
if ~iscell(I)
69+
70+
% extracts initial and final times
71+
t0 = I(1);
72+
tf = I(2);
73+
74+
% makes step size negative if t0 > tf
75+
if (t0 > tf)
76+
h = -h;
77+
end
78+
79+
% defines condition function
80+
C = @(t,y) t <= tf;
81+
82+
% indicates that final time is known
83+
final_time_known = true;
84+
85+
% number of subintervals between sample times
86+
N = ceil((tf-t0)/h);
87+
88+
% turns waitbar on with custom message
89+
if (nargin == 6) && ischar(wb)
90+
display_waitbar = true;
91+
[wb,prop] = initialize_waitbar(msg);
92+
93+
% turns waitbar on with default message
94+
elseif (nargin == 6) && wb
95+
display_waitbar = true;
96+
[wb,prop] = initialize_waitbar('Solving ODE...');
97+
98+
% turns waitbar off otherwise
99+
else
100+
display_waitbar = false;
101+
102+
end
103+
104+
% ---------------------
105+
% Event detection mode.
106+
% ---------------------
107+
108+
else
109+
110+
% extracts initial time and condition function
111+
t0 = I{1};
112+
C = I{2};
113+
114+
% indicates that final time is unknown
115+
final_time_known = false;
116+
117+
end
118+
119+
% -------------------
120+
% Integration method.
121+
% -------------------
122+
123+
% defaults integration method to RK4
124+
if (nargin < 5) || isempty(method)
125+
method = 'RK4';
126+
end
127+
128+
% sets propagation function for single-step methods
129+
if strcmpi(method,'RK1_euler')
130+
propagate = @(t,y) RK1_euler(f,t,y,h);
131+
elseif strcmpi(method,'RK2')
132+
propagate = @(t,y) RK2(f,t,y,h);
133+
elseif strcmpi(method,'RK2_heun')
134+
propagate = @(t,y) RK2_heun(f,t,y,h);
135+
elseif strcmpi(method,'RK2_ralston')
136+
propagate = @(t,y) RK2_ralston(f,t,y,h);
137+
elseif strcmpi(method,'RK3')
138+
propagate = @(t,y) RK3(f,t,y,h);
139+
elseif strcmpi(method,'RK3_heun')
140+
propagate = @(t,y) RK3_heun(f,t,y,h);
141+
elseif strcmpi(method,'RK3_ralston')
142+
propagate = @(t,y) RK3_ralston(f,t,y,h);
143+
elseif strcmpi(method,'SSPRK3')
144+
propagate = @(t,y) SSPRK3(f,t,y,h);
145+
elseif strcmpi(method,'RK4')
146+
propagate = @(t,y) RK4(f,t,y,h);
147+
elseif strcmpi(method,'RK4_38')
148+
propagate = @(t,y) RK4_38(f,t,y,h);
149+
elseif strcmpi(method,'RK4_ralston')
150+
propagate = @(t,y) RK4_ralston(f,t,y,h);
151+
end
152+
153+
% sets propagation function for multi-step predictor methods
154+
if strcmpi(method,'AB2')
155+
propagate = @(t,F) AB2(f,t,F,h);
156+
elseif strcmpi(method,'AB3')
157+
propagate = @(t,F) AB3(f,t,F,h);
158+
elseif strcmpi(method,'AB4')
159+
propagate = @(t,F) AB5(f,t,F,h);
160+
elseif strcmpi(method,'AB3')
161+
propagate = @(t,F) AB3(f,t,F,h);
162+
elseif strcmpi(method,'AB3')
163+
propagate = @(t,F) AB3(f,t,F,h);
164+
elseif strcmpi(method,'AB3')
165+
propagate = @(t,F) AB3(f,t,F,h);
166+
elseif strcmpi(method,'AB3')
167+
propagate = @(t,F) AB3(f,t,F,h);
168+
elseif strcmpi(method,'AB3')
169+
propagate = @(t,F) AB3(f,t,F,h);
170+
end
171+
172+
% determines if the method is a single-step method
173+
if strcmpi(method(1:2),'RK')
174+
single_step = true;
175+
else
176+
single_step = false;
177+
end
178+
179+
% determines order for multistep methods
180+
if ~single_step
181+
m = str2double(method(end));
182+
end
183+
184+
% --------------------------------------------------
185+
% Preallocates arrays and stores initial conditions.
186+
% --------------------------------------------------
187+
188+
% preallocates time vector and solution matrix
189+
t = zeros(10000,1);
190+
y = zeros(length(y0),length(t));
191+
192+
% stores initial conditions
193+
t(1) = t0;
194+
y(:,1) = y0;
195+
196+
% ---------------------------------------
197+
% Solves ODE (using single-step methods).
198+
% ---------------------------------------
199+
200+
if single_step
201+
202+
% state vector propagation while condition is satisfied
203+
n = 1;
204+
while C(t(n),y(:,n))
205+
206+
% expands t and y if needed
207+
if (n+1) > length(t)
208+
[t,y] = expand_solution_arrays(t,y);
209+
end
210+
211+
% state vector propagated to next sample time
212+
y(:,n+1) = propagate(t(n),y(:,n));
213+
214+
% increments time and loop index
215+
t(n+1) = t(n)+h;
216+
n = n+1;
217+
218+
% updates waitbar
219+
if final_time_known && display_waitbar
220+
prop = update_waitbar(n,N,wb,prop);
221+
end
222+
223+
end
224+
225+
% --------------------------------------
226+
% Solves ODE (using multi-step methods).
227+
% --------------------------------------
228+
229+
else
230+
231+
% time vector for first m sample times
232+
t(1:m) = (t0:h:(t0+(m-1)*h))';
233+
234+
% propagating state vector using RK4 for first "m" sample times
235+
for n = 1:m
236+
y(:,n+1) = RK4(f,t(n),y(:,n),h);
237+
end
238+
239+
% initializes F matrix (stores function evaluations for first m
240+
% sample times in first m columns, state vector at (m+1)th sample
241+
% time in (m+1)th column)
242+
F = zeros(length(y0),m+1);
243+
for n = 1:m
244+
F(:,n) = f(t(n),y(:,n));
245+
end
246+
F(:,n+1) = y(:,m+1);
247+
248+
% state vector propagation while condition is satisfied
249+
n = 1;
250+
while C(t(n),y(:,n))
251+
252+
% expands t and y if needed
253+
if (n+1) > length(t)
254+
[t,y] = expand_solution_arrays(t,y);
255+
end
256+
257+
% updates F matrix
258+
F = propagate(t(n),F);
259+
260+
% extracts/stores state vector propagated to next sample time
261+
y(:,n+1) = F(:,m+1);
262+
263+
% increments time and loop index
264+
t(n+1) = t(n)+h;
265+
n = n+1;
266+
267+
end
268+
269+
end
270+
271+
% -----------------
272+
% Final formatting.
273+
% -----------------
274+
275+
% trims arrays
276+
y = y(:,1:(n-1));
277+
t = t(1:(n-1));
278+
279+
% transposes solution matrix so it is returned in "standard form"
280+
y = y';
281+
282+
% closes waitbar
283+
if final_time_known && display_waitbar
284+
close(wb);
285+
end
286+
287+
end

0 commit comments

Comments
 (0)