Skip to content

Commit 737b840

Browse files
tamaskistamaskis
committed
fixed typos
1 parent 4e30143 commit 737b840

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

46 files changed

+9005
-0
lines changed

‎toolbox/adamsbashforth/AB2.m‎

Lines changed: 317 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,317 @@
1+
%==========================================================================
2+
%
3+
% AB2 Adams-Bashforth 2nd-order method.
4+
%
5+
% [t,y] = AB2(f,[t0,tf],y0,h)
6+
% [t,y] = AB2(f,{t0,C},y0,h)
7+
% [t,y] = AB2(__,wb)
8+
%
9+
% See also AB3, AB4, AB5, AB6, AB7, AB8.
10+
%
11+
% Copyright © 2021 Tamas Kis
12+
% Last Update: 2021-12-22
13+
% Website: https://tamaskis.github.io
14+
% Contact: tamas.a.kis@outlook.com
15+
%
16+
% TECHNICAL DOCUMENTATION:
17+
% https://tamaskis.github.io/documentation/Fixed_Step_ODE_Solvers.pdf
18+
%
19+
% REFERENCES:
20+
% [1] https://en.wikipedia.org/wiki/Linear_multistep_method
21+
% [2] Montenbruck and Gill, "Satellite Orbits" (pp. 135-136)
22+
%
23+
%--------------------------------------------------------------------------
24+
%
25+
% ------
26+
% INPUT:
27+
% ------
28+
% f - (1×1 function_handle) dy/dt = f(t,y) --> multivariate,
29+
% vector-valued function (f:R×Rp->Rp) 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, t0, and function
33+
% handle for condition function, C(t,y)
34+
% (C:R×Rp->B)
35+
% y0 - (p×1 double) initial condition
36+
% h - (1×1 double) step size
37+
% wb - (OPTIONAL) (1×1 logical or char) waitbar parameters
38+
% --> input as "true" if you want waitbar with default
39+
% message displayed
40+
% --> input as a char array storing a message if you want a
41+
% custom message displayed on the waitbar
42+
%
43+
% -------
44+
% OUTPUT:
45+
% -------
46+
% t - ((N+1)×1 double) time vector
47+
% y - ((N+1)×p double) matrix storing time history of state vector
48+
%
49+
% -----
50+
% NOTE:
51+
% -----
52+
% --> p = dimension of state vector (for the scalar case, p = 1)
53+
% --> m = order of multistep method
54+
% --> N+1 = length of time vector
55+
% --> The ith row of "y" is the TRANSPOSE of the state vector (i.e. the
56+
% solution corresponding to the ith time in "t"). This convention is
57+
% chosen to match the convention used by MATLAB's ODE suite.
58+
%
59+
%==========================================================================
60+
function [t,y] = AB2(f,I,y0,h,wb)
61+
62+
% ------------------------------
63+
% Order of the multistep method.
64+
% ------------------------------
65+
66+
m = 2;
67+
68+
% -------------------
69+
% Setting up waitbar.
70+
% -------------------
71+
72+
% determines if waitbar is on or off
73+
if (nargin < 5) || (islogical(wb) && ~wb)
74+
display_waitbar = false;
75+
else
76+
display_waitbar = true;
77+
end
78+
79+
% sets the waitbar message (defaults to 'Solving ODE...')
80+
if display_waitbar
81+
if ischar(wb)
82+
msg = wb;
83+
else
84+
msg = 'Solving ODE...';
85+
end
86+
end
87+
88+
% initialize cutoff proportion needed to trigger waitbar update to 0.1
89+
if display_waitbar, prop = 0.1; end
90+
91+
% ---------------------------------------
92+
% Determines which implementation to use.
93+
% ---------------------------------------
94+
95+
if iscell(I)
96+
t0 = I{1};
97+
C = I{2};
98+
time_detection_implementation = false;
99+
else
100+
t0 = I(1);
101+
tf = I(2);
102+
time_detection_implementation = true;
103+
end
104+
105+
% --------------------------------------------------------
106+
% Time detection implementation (solves until final time).
107+
% --------------------------------------------------------
108+
109+
if time_detection_implementation
110+
111+
% initializes the waitbar
112+
if display_waitbar, wb = waitbar(0,msg); end
113+
114+
% makes step size negative if t0 > tf
115+
if t0 > tf
116+
h = -h;
117+
end
118+
119+
% number of subintervals between sample times
120+
N = ceil((tf-t0)/h);
121+
122+
% last element of the time vector
123+
tN = t0+N*h;
124+
125+
% defines time vector, preallocates solution matrix, and
126+
% preallocates matrix to store previous m function evaluations
127+
t = (t0:h:tN)';
128+
y = zeros(length(y0),length(t));
129+
fm = zeros(length(y0),m);
130+
131+
% stores initial condition in solution matrix
132+
y(:,1) = y0;
133+
134+
% propagating state vector using RK4 for first "m" sample times
135+
for n = 1:m
136+
137+
% current sample time and state vector
138+
tn = t(n);
139+
yn = y(:,n);
140+
141+
% k terms
142+
k1 = f(tn,yn);
143+
k2 = f(tn+h/2,yn+h*k1/2);
144+
k3 = f(tn+h/2,yn+h*k2/2);
145+
k4 = f(tn+h,yn+h*k3);
146+
147+
% state vector propagated to next sample time
148+
y(:,n+1) = yn+(h/6)*(k1+2*k2+2*k3+k4);
149+
150+
% updates waitbar
151+
if display_waitbar, prop = update_waitbar(n,N,wb,prop); end
152+
153+
end
154+
155+
% stores function evaluations for first m sample times
156+
for n = 1:m
157+
fm(:,n) = f(t(n),y(:,n));
158+
end
159+
160+
% propagating state vector using AB2
161+
for n = (m+1):N
162+
163+
% updates stored function evaluations
164+
fm = [fm(:,2:end),f(t(n),y(:,n))];
165+
166+
% state vector propagated to next sample time
167+
y(:,n+1) = y(:,n)+(h/2)*(3*fm(:,m)-fm(:,m-1));
168+
169+
% updates waitbar
170+
if display_waitbar, prop = update_waitbar(n,N,wb,prop); end
171+
172+
end
173+
174+
% linearly interpolates for solution at tf
175+
y(:,N+1) = y(:,N)+((y(:,N+1)-y(:,N))/(t(N+1)-t(N)))*(tf-t(N));
176+
177+
% replaces last element of "t" with tf
178+
t(N+1) = tf;
179+
180+
% closes waitbar
181+
if display_waitbar, close(wb); end
182+
183+
% ---------------------------------------------------------------------
184+
% Event detection implementation (solves while condition is satisfied).
185+
% ---------------------------------------------------------------------
186+
187+
else
188+
189+
% preallocates time vector, solution matrix, and array to store
190+
% previous m function evaluations
191+
t = zeros(10000,1);
192+
y = zeros(length(y0),length(t));
193+
fm = zeros(length(y0),m);
194+
195+
% time vector for first m sample times
196+
t(1:m) = (t0:h:(t0+(m-1)*h))';
197+
198+
% stores initial condition in solution matrix
199+
t(1) = t0;
200+
y(:,1) = y0;
201+
202+
% propagating state vector using RK4 for first m sample times
203+
for n = 1:m
204+
205+
% current sample time and state vector
206+
tn = t(n);
207+
yn = y(:,n);
208+
209+
% k terms
210+
k1 = f(tn,yn);
211+
k2 = f(tn+h/2,yn+k1*h/2);
212+
k3 = f(tn+h/2,yn+k2*h/2);
213+
k4 = f(tn+h,yn+k3*h);
214+
215+
% state vector propagated to next sample time
216+
y(:,n+1) = yn+(h/6)*(k1+2*k2+2*k3+k4);
217+
218+
end
219+
220+
% stores function evaluations for first m sample times
221+
for n = 1:m
222+
fm(:,n) = f(t(n),y(:,n));
223+
end
224+
225+
% state vector propagation while condition is satisfied
226+
n = m;
227+
while C(t(n),y(:,n))
228+
229+
% expands t and y if needed
230+
if (n+1) > length(t)
231+
[t,y] = expand_solution_arrays(t,y);
232+
end
233+
234+
% updates stored function evaluations
235+
fm = [fm(:,2:end),f(t(n),y(:,n))];
236+
237+
% state vector propagated to next sample time
238+
y(:,n+1) = y(:,n)+(h/2)*(3*fm(:,m)-fm(:,m-1));
239+
240+
% increments time and loop index
241+
t(n+1) = t(n)+h;
242+
n = n+1;
243+
244+
end
245+
246+
% trims arrays
247+
y = y(:,1:(n-1));
248+
t = t(1:(n-1));
249+
250+
end
251+
252+
% transposes solution matrix so it is returned in "standard form"
253+
y = y';
254+
255+
% -------------
256+
% Subfunctions.
257+
% -------------
258+
259+
%----------------------------------------------------------------------
260+
% expand_solution_arrays Expands the arrays storing the ODE solution.
261+
% This function is used to preallocate more space for the solution
262+
% arrays (i.e. the time vector and the solution matrix) once they have
263+
% become full.
264+
%----------------------------------------------------------------------
265+
%
266+
% INPUT:
267+
% t - ((N+1)×1 double) time vector
268+
% y - (p×(N+1) double) solution matrix
269+
%
270+
% OUTPUT:
271+
% t_new - (2(N+1)×1 double) expanded time vector
272+
% y_new - (p×2(N+1) double) expanded solution matrix
273+
%
274+
% NOTE:
275+
% --> p = dimension of state vector (for the scalar case, p = 1)
276+
% --> N+1 = original length of time vector
277+
%
278+
%----------------------------------------------------------------------
279+
function [t_new,y_new] = expand_solution_arrays(t,y)
280+
t_new = [t;zeros(length(t),1)];
281+
y_new = [y,zeros(size(y,1),length(t))];
282+
end
283+
284+
%----------------------------------------------------------------------
285+
% update_waitbar Updates the waitbar.
286+
%----------------------------------------------------------------------
287+
%
288+
% INPUT:
289+
% wb - (1×1 Figure) waitbar
290+
% n - (1×1 double) current sample number (i.e. iteration)
291+
% N - (1×1 double) total number of samples (i.e. iterations)
292+
% prop - (1×1 double) cutoff proportion to trigger waitbar update
293+
%
294+
% OUTPUT:
295+
% prop - (1×1 double) cutoff proportion to trigger waitbar update
296+
%
297+
% NOTE:
298+
% --> "prop" is an integer multiple of 0.1 so that the waitbar is
299+
% only updated after every additional 10% of progress.
300+
%
301+
%----------------------------------------------------------------------
302+
function prop = update_waitbar(i,N,wb,prop)
303+
304+
% only updates waitbar if current proportion exceeds cutoff prop.
305+
if i/N > prop
306+
307+
% updates waitbar
308+
waitbar(i/N,wb);
309+
310+
% updates cutoff proportion needed to trigger waitbar update
311+
prop = prop+0.1;
312+
313+
end
314+
315+
end
316+
317+
end

0 commit comments

Comments
 (0)