Mercurial > repos > public > sbplib
annotate +time/Rungekutta4SecondOrder.m @ 1341:663eb90a4559 feature/D2_boundary_opt
Pass logic grid along to diracDiscr
author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
---|---|
date | Fri, 22 Jul 2022 16:37:49 +0200 |
parents | 74eec7e69b63 |
children |
rev | line source |
---|---|
0 | 1 classdef Rungekutta4SecondOrder < time.Timestepper |
2 properties | |
3 F | |
4 k | |
5 t | |
6 w | |
7 m | |
8 D | |
9 E | |
10 S | |
11 M | |
12 C | |
13 n | |
14 end | |
15 | |
16 | |
17 methods | |
549
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
18 % Solves u_tt = Du + Eu_t + S by |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
19 % Rewriting on first order form: |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
20 % w_t = M*w + C(t) |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
21 % where |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
22 % M = [ |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
23 % 0, I; |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
24 % D, E; |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
25 % ] |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
26 % and |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
27 % C(t) = [ |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
28 % 0; |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
29 % S(t) |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
30 % ] |
222
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
31 % D, E, S can either all be constants or all be function handles, |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
32 % They can also be omitted by setting them equal to the empty matrix. |
0 | 33 function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t) |
34 obj.D = D; | |
35 obj.E = E; | |
36 obj.S = S; | |
37 obj.m = length(v0); | |
13
b18d3d201a71
Fixed initialization of step counter in timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
38 obj.n = 0; |
1258
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
39 |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
40 % Construct RHS system on first order form. |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
41 % 3 different setups are handeled: |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
42 % If none of D, E, S are time-dependent (e.g function_handles) |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
43 % then a complete matrix of the RHS is constructed. |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
44 % If the source term S is time-depdentent, then the first order system matrix M |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
45 % is formed via D and E, while S is kept as a function handle |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
46 % If D or E are time-dependent, then all terns are treated as function handles. |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
47 if ~isa(D, 'function_handle') && ~isa(E, 'function_handle') && isa(S, 'function_handle') |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
48 default_arg('D', sparse(obj.m, obj.m)); |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
49 default_arg('E', sparse(obj.m, obj.m)); |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
50 default_arg('S', @(t)sparse(obj.m, 1)); |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
51 |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
52 I = speye(obj.m); |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
53 O = sparse(obj.m,obj.m); |
222
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
54 |
1258
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
55 obj.M = [ |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
56 O, I; |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
57 D, E; |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
58 ]; |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
59 |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
60 obj.k = k; |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
61 obj.t = t0; |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
62 obj.w = [v0; v0t]; |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
63 obj.F = @(w,t)obj.rhsTimeDepS(w,t,obj.M,S,obj.m); |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
64 elseif isa(D, 'function_handle') || isa(E, 'function_handle') |
222
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
65 default_arg('D', @(t)sparse(obj.m, obj.m)); |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
66 default_arg('E', @(t)sparse(obj.m, obj.m)); |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
67 default_arg('S', @(t)sparse(obj.m, 1) ); |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
68 |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
69 if ~isa(D, 'function_handle') |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
70 D = @(t)D; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
71 end |
345
7b5ef8b89268
Rungekutta bug fixes.
Jonatan Werpers <jonatan@werpers.com>
parents:
222
diff
changeset
|
72 if ~isa(E, 'function_handle') |
222
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
73 E = @(t)E; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
74 end |
345
7b5ef8b89268
Rungekutta bug fixes.
Jonatan Werpers <jonatan@werpers.com>
parents:
222
diff
changeset
|
75 if ~isa(S, 'function_handle') |
222
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
76 S = @(t)S; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
77 end |
0 | 78 |
222
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
79 obj.k = k; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
80 obj.t = t0; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
81 obj.w = [v0; v0t]; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
82 |
549
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
83 % Avoid matrix formulation because it is VERY slow |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
84 obj.F = @(w,t)[ |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
85 w(obj.m+1:end); |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
86 D(t)*w(1:obj.m) + E(t)*w(obj.m+1:end) + S(t); |
ae905a11e32c
Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents:
345
diff
changeset
|
87 ]; |
0 | 88 else |
222
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
89 default_arg('D', sparse(obj.m, obj.m)); |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
90 default_arg('E', sparse(obj.m, obj.m)); |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
91 default_arg('S', sparse(obj.m, 1) ); |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
92 |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
93 I = speye(obj.m); |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
94 O = sparse(obj.m,obj.m); |
0 | 95 |
222
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
96 obj.M = [ |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
97 O, I; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
98 D, E; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
99 ]; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
100 obj.C = [ |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
101 zeros(obj.m,1); |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
102 S; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
103 ]; |
0 | 104 |
222
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
105 obj.k = k; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
106 obj.t = t0; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
107 obj.w = [v0; v0t]; |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
108 |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
109 obj.F = @(w,t)(obj.M*w + obj.C); |
e7e73173d44d
time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
110 end |
0 | 111 end |
112 | |
113 function [v,t] = getV(obj) | |
114 v = obj.w(1:end/2); | |
115 t = obj.t; | |
116 end | |
117 | |
118 function [vt,t] = getVt(obj) | |
119 vt = obj.w(end/2+1:end); | |
120 t = obj.t; | |
121 end | |
122 | |
123 function obj = step(obj) | |
124 obj.w = time.rk4.rungekutta_4(obj.w, obj.t, obj.k, obj.F); | |
125 obj.t = obj.t + obj.k; | |
126 obj.n = obj.n + 1; | |
127 end | |
128 end | |
129 | |
130 | |
131 methods (Static) | |
132 function k = getTimeStep(lambda) | |
133 k = rk4.get_rk4_time_step(lambda); | |
134 end | |
1258
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
135 |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
136 function F = rhsTimeDepS(w,t,M,S,m) |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
137 F = M*w; |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
138 F(m+1:end) = F(m+1:end) + S(t); |
74eec7e69b63
Speed up computations using RungeKuttaSecondOrder when only the source term is time dependent.
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
549
diff
changeset
|
139 end |
0 | 140 end |
141 | |
142 end |