Mercurial > repos > public > sbplib

annotate +time/Rungekutta4SecondOrder.m @ 887:50d5a3843099 feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Rename package rk4 to rk
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 15 Nov 2018 16:42:58 -0800
parents 8894e9c49e40
children f5e14e5986b5
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
1 classdef Rungekutta4SecondOrder < time.Timestepper
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2 properties
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
3 F
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
4 k
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
5 t
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6 w
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
7 m
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8 D
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9 E
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
10 S
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
11 M
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
12 C
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
13 n
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
14 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
15
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
16
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
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
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
33 function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
34 obj.D = D;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
35 obj.E = E;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
36 obj.S = S;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
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;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
39
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
40
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
41 if isa(D, 'function_handle') || isa(E, 'function_handle') || isa(S, '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
42 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
43 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
44 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
45
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
46 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
47 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
48 end
345
7b5ef8b89268 Rungekutta bug fixes.
Jonatan Werpers <jonatan@werpers.com>
parents: 222
diff changeset
49 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
50 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
51 end
345
7b5ef8b89268 Rungekutta bug fixes.
Jonatan Werpers <jonatan@werpers.com>
parents: 222
diff changeset
52 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
53 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
54 end
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
55
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
56 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
57 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
58 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
59
549
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
60 % 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
61 obj.F = @(w,t)[
ae905a11e32c Change away from matrix formulation in Rk4 second order form
Jonatan Werpers <jonatan@werpers.com>
parents: 345
diff changeset
62 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
63 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
64 ];
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
65 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
66
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('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
68 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
69 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
70
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 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
72 O = sparse(obj.m,obj.m);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
73
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
74 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
75 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
76 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
77 ];
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
78 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
79 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
80 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
81 ];
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
82
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
83 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
84 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
85 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
86
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
87 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
88 end
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
89 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
90
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
91 function [v,t] = getV(obj)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
92 v = obj.w(1:end/2);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
93 t = obj.t;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
94 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
95
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
96 function [vt,t] = getVt(obj)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
97 vt = obj.w(end/2+1:end);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
98 t = obj.t;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
99 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
100
816
b5e5b195da1e Add getState to timesteppers, returning the relevant state of the timestepper
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 13
diff changeset
101 function state = getState(obj)
b5e5b195da1e Add getState to timesteppers, returning the relevant state of the timestepper
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 13
diff changeset
102 [v, t] = obj.getV();
b5e5b195da1e Add getState to timesteppers, returning the relevant state of the timestepper
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 13
diff changeset
103 [vt] = obj.getVt();
b5e5b195da1e Add getState to timesteppers, returning the relevant state of the timestepper
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 13
diff changeset
104 state = struct('v', v, 'vt', vt, 't', t, 'k', obj.k);
b5e5b195da1e Add getState to timesteppers, returning the relevant state of the timestepper
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 13
diff changeset
105 end
b5e5b195da1e Add getState to timesteppers, returning the relevant state of the timestepper
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 13
diff changeset
106
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
107 function obj = step(obj)
887
50d5a3843099 Rename package rk4 to rk
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents: 886
diff changeset
108 obj.w = time.rk.rungekutta_4(obj.w, obj.t, obj.k, obj.F);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
109 obj.t = obj.t + obj.k;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
110 obj.n = obj.n + 1;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
111 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
112 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
113
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
114
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
115 methods (Static)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
116 function k = getTimeStep(lambda)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
117 k = rk4.get_rk4_time_step(lambda);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
118 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
119 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
120
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
121 end