Skip to content

Commit 012fe0b

Browse files
committed
updated examples, reorganized functions, updated solve_ivp
1 parent d091291 commit 012fe0b

File tree

9 files changed

+181
-74
lines changed

9 files changed

+181
-74
lines changed

examples/example_brute_force.m

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@
44
% Example of a "brute force" method for solving IVPs.
55
%
66
% Copyright © 2021 Tamas Kis
7-
% Last Update: 2022-06-05
7+
% Last Update: 2022-06-06
88
% Website: https://tamaskis.github.io
99
% Contact: tamas.a.kis@outlook.com
1010

@@ -62,7 +62,7 @@
6262
n = n+1;
6363

6464
end
65-
65+
6666
% trims arrays
6767
t = t(1:(n-1));
6868
x = x(1:(n-1));

examples/example_define_ivp.m

Lines changed: 75 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,75 @@
1+
%% example_define_ivp.m
2+
% IVP Solver Toolbox
3+
%
4+
% Example of defining an IVP.
5+
%
6+
% Copyright © 2021 Tamas Kis
7+
% Last Update: 2022-06-06
8+
% Website: https://tamaskis.github.io
9+
% Contact: tamas.a.kis@outlook.com
10+
11+
12+
13+
%% SCRIPT SETUP
14+
15+
% clears Workspace and Command Window, closes all figures
16+
clear; clc; close all;
17+
18+
19+
20+
%% SOLUTION
21+
22+
% parameters
23+
b = 5; % damping constant [N.s/m]
24+
k = 1; % spring constant [N/m]
25+
m = 2; % mass [kg]
26+
x0 = 1; % initial position [m]
27+
dx0 = 0; % initial velocity [m/s]
28+
29+
% forcing function
30+
F = @(t) cos(pi*t);
31+
32+
% initial condition
33+
y0 = [x0;
34+
dx0];
35+
36+
% ODE
37+
f = @(t,y) [y(2);
38+
-(b/m)*y(2)-(k/m)*y(1)+(1/m)*F(t)];
39+
40+
41+
42+
%% ALTERNATE SOLUTION
43+
44+
% parameters
45+
b = 5; % damping constant [N.s/m]
46+
k = 1; % spring constant [N/m]
47+
m = 2; % mass [kg]
48+
x0 = 1; % initial position [m]
49+
dx0 = 0; % initial velocity [m/s]
50+
51+
% forcing function
52+
F = @(t) cos(pi*t);
53+
54+
% initial condition
55+
y0 = [x0;
56+
dx0];
57+
58+
% assigns function handle to ODE
59+
f = @(t,y) f_extra(t,y,b,k,m,F);
60+
61+
% defines ODE
62+
function dy = f_extra(t,y,b,k,m,F)
63+
64+
% unpacks state vector
65+
x = y(1);
66+
xdot = y(2);
67+
68+
% preallocates state vector derivative
69+
dy = zeros(size(y));
70+
71+
% assembles state vector derivative
72+
dy(1) = xdot;
73+
dy(2) = -(b/m)*xdot-(k/m)*x+(1/m)*F(t);
74+
75+
end

examples/example_lorenz.m

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@
55
% the thumbnail image for the IVP solver toolbox.
66
%
77
% Copyright © 2021 Tamas Kis
8-
% Last Update: 2022-06-05
8+
% Last Update: 2022-06-06
99
% Website: https://tamaskis.github.io
1010
% Contact: tamas.a.kis@outlook.com
1111

@@ -31,7 +31,7 @@
3131
x(1)*x(2)-beta*x(3)];
3232

3333
% solution
34-
[t,y] = odeabm(f,[0,100],[0;1;1.05],0.001);
34+
[t,y] = solve_ivp(f,[0,100],[0;1;1.05],0.001,'ABM8');
3535

3636
% plots figure
3737
figure;

examples/example_matrix_ode.m

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@
55
% equation).
66
%
77
% Copyright © 2021 Tamas Kis
8-
% Last Update: 2022-06-05
8+
% Last Update: 2022-06-06
99
% Website: https://tamaskis.github.io
1010
% Contact: tamas.a.kis@outlook.com
1111

@@ -57,14 +57,14 @@
5757
f = odefun_mat2vec(F);
5858

5959
% final condition
60-
yT = odeIC_mat2vec(PT);
60+
yT = ivpIC_mat2vec(PT);
6161

6262
% solves vector-valued IVP using a step size of h = 0.001
63-
[~,y] = RK4(f,[T,0],yT,0.001);
63+
[~,y] = solve_ivp(f,[T,0],yT,0.001,'RK4');
6464

6565
% transforms solution matrix for vector-valued IVP into solution array for
6666
% matrix-valued IVP
67-
P = odesol_vec2mat(y);
67+
P = ivpsol_vec2mat(y);
6868

6969
% solution for P0 (will be at end of array since P solved for backwards in
7070
% time)
@@ -84,9 +84,9 @@
8484
% store initial condition
8585
P(:,:,1) = PT;
8686

87-
% solving using "RK4_step"
87+
% solving using "RK4"
8888
for i = 1:(length(t)-1)
89-
P(:,:,i+1) = RK4_step(F,t(i),P(:,:,i),h);
89+
P(:,:,i+1) = RK4(F,t(i),P(:,:,i),h);
9090
end
9191

9292
% solution for P0 using one-step propagation

examples/example_pass_extra.m

Lines changed: 51 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,10 @@
11
%% example_pass_extra.m
22
% IVP Solver Toolbox
33
%
4-
% Example of passing extra parameters to a function defining an IVP.
4+
% Example of passing extra parameters to functions.
55
%
66
% Copyright © 2021 Tamas Kis
7-
% Last Update: 2022-06-05
7+
% Last Update: 2022-06-06
88
% Website: https://tamaskis.github.io
99
% Contact: tamas.a.kis@outlook.com
1010

@@ -17,59 +17,65 @@
1717

1818

1919

20-
%% SOLUTION
20+
%% EXAMPLE #1
2121

22-
% parameters
23-
b = 5; % damping constant [N.s/m]
24-
k = 1; % spring constant [N/m]
25-
m = 2; % mass [kg]
26-
x0 = 1; % initial position [m]
27-
dx0 = 0; % initial velocity [m/s]
22+
% function definition
23+
a = 5; b = 5; c = 5; d = 5;
24+
f = @(x,y) -a*(x-b)^2-c*(y-d)^2;
2825

29-
% forcing function
30-
F = @(t) cos(pi*t);
26+
% call function, update "a", then call function again
27+
f(2,2)
28+
a = 20;
29+
f(2,2)
3130

32-
% initial condition
33-
y0 = [x0;
34-
dx0];
31+
% NOTE: Both evaluations of "f" yield the same result.
3532

36-
% differential equation
37-
f = @(t,y) [y(2);
38-
-(b/m)*y(2)-(k/m)*y(1)+(1/m)*F(t)];
3933

4034

35+
%% EXAMPLE #2
4136

42-
%% ALTERNATE SOLUTION
37+
% original function definition
38+
a = 5; b = 5; c = 5; d = 5;
39+
f = @(x,y) -a*(x-b)^2-c*(y-d)^2;
40+
f(2,2)
4341

44-
% parameters
45-
b = 5; % damping constant [N.s/m]
46-
k = 1; % spring constant [N/m]
47-
m = 2; % mass [kg]
48-
x0 = 1; % initial position [m]
49-
dx0 = 0; % initial velocity [m/s]
42+
% update values of constants
43+
a = 10; b = 10; c = 10; d = 10;
44+
f = @(x,y) -a*(x-b)^2-c*(y-d)^2;
45+
f(2,2)
5046

51-
% forcing function
52-
F = @(t) cos(pi*t);
47+
% NOTE: The two evaluations of "f" no longer yield the same result.
5348

54-
% initial condition
55-
y0 = [x0;
56-
dx0];
5749

58-
% assigns function handle to differential equation
59-
f = @(t,y) f_extra(t,y,b,k,m,F);
6050

61-
% defines differential equation
62-
function dy = f_extra(t,y,b,k,m,F)
51+
%% EXAMPLE #3
6352

64-
% unpacks state vector
65-
x = y(1);
66-
xdot = y(2);
67-
68-
% preallocates state vector derivative
69-
dy = zeros(size(y));
70-
71-
% assembles state vector derivative
72-
dy(1) = xdot;
73-
dy(2) = -(b/m)*xdot-(k/m)*x+(1/m)*F(t);
74-
75-
end
53+
% define function where constants can vary as well
54+
f_extra = @(x,y,a,b,c,d) -a*(x-b)^2-c*(y-d)^2;
55+
56+
% define original function by assigning function handle to f_extra
57+
f = @(x,y) f_extra(x,y,5,5,5,5);
58+
f(2,2)
59+
60+
% update values of constants
61+
f = @(x,y) f_extra(x,y,10,10,10,10);
62+
f(2,2)
63+
64+
% NOTE: This example is an alternate solution to Example #2.
65+
66+
67+
68+
%% EXAMPLE #4
69+
70+
% sets values of constants
71+
a = 10; b = 10; c = 10; d = 10;
72+
f = @(x,y) f_extra2(x,y,a,b,c,d);
73+
f(2,2)
74+
75+
% MATLAB function must be declared at end
76+
function f = f_extra2(x,y,a,b,c,d)
77+
f = -a*(x-b)^2-c*(y-d)^2;
78+
end
79+
80+
% NOTE: This example is an alternate solution to the second part of
81+
% Example #2.

examples/example_return_extra.m

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@
44
% Example of returning extra parameters from an IVP solution.
55
%
66
% Copyright © 2021 Tamas Kis
7-
% Last Update: 2022-06-05
7+
% Last Update: 2022-06-06
88
% Website: https://tamaskis.github.io
99
% Contact: tamas.a.kis@outlook.com
1010

examples/example_waitbar.m

Lines changed: 6 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@
44
% Examples demonstrating how to use waitbar.
55
%
66
% Copyright © 2021 Tamas Kis
7-
% Last Update: 2022-06-05
7+
% Last Update: 2022-06-06
88
% Website: https://tamaskis.github.io
99
% Contact: tamas.a.kis@outlook.com
1010

@@ -33,12 +33,13 @@
3333

3434
%% EXAMPLE #1 - WAITBAR WITH DEFAULT MESSAGE
3535

36-
% solves IVP, displaying waitbar with default message
37-
[t,y] = ABM8(f,[0,200],[0;1;1.05],0.001,true);
36+
% solves IVP using RK4, displaying waitbar with default message
37+
[t,y] = solve_ivp(f,[0,200],[0;1;1.05],0.001,[],true);
3838

3939

4040

4141
%% EXAMPLE #2 - WAITBAR WITH CUSTOM MESSAGE
4242

43-
% solves IVP, displaying waitbar with custom message
44-
[t,y] = ABM8(f,[0,200],[0;1;1.05],0.001,'Solving Lorenz system...');
43+
% solves IVP using ABM8, displaying waitbar with custom message
44+
[t,y] = solve_ivp(f,[0,200],[0;1;1.05],0.001,'ABM8',...
45+
'Solving Lorenz system...');
Lines changed: 14 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@
55
% [t_new,y_new] = expand_solution_arrays(t,y)
66
%
77
% Copyright © 2021 Tamas Kis
8-
% Last Update: 2022-06-05
8+
% Last Update: 2022-06-06
99
% Website: https://tamaskis.github.io
1010
% Contact: tamas.a.kis@outlook.com
1111
%
@@ -40,6 +40,17 @@
4040
%
4141
%==========================================================================
4242
function [t_new,y_new] = expand_solution_arrays(t,y)
43-
t_new = [t;zeros(length(t),1)];
44-
y_new = [y,zeros(size(y,1),length(t))];
43+
44+
% number of subintervals
45+
N = length(t)-1;
46+
47+
% state dimension
48+
p = size(y,1);
49+
50+
% expands time vector
51+
t_new = [t;zeros(N+1,1)];
52+
53+
% expands solution matrix
54+
y_new = [y,zeros(p,N+1)];
55+
4556
end

0 commit comments

Comments
 (0)