Mercurial > repos > public > sbplib

annotate +time/Rungekutta4SecondOrder.m @ 345:7b5ef8b89268 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Rungekutta bug fixes.
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 16 Nov 2016 15:39:16 -0800
parents e7e73173d44d
children ae905a11e32c
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
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
18 % Solves u_tt = Du + Eu_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
19 % 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
20 % 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
21 function obj = Rungekutta4SecondOrder(D, E, S, k, t0, v0, v0t)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
22 obj.D = D;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
23 obj.E = E;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
24 obj.S = S;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
25 obj.m = length(v0);
13
b18d3d201a71 Fixed initialization of step counter in timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents: 0
diff changeset
26 obj.n = 0;
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
27
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
28
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
29 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
30 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
31 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
32 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
33
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
34 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
35 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
36 end
345
7b5ef8b89268 Rungekutta bug fixes.
Jonatan Werpers <jonatan@werpers.com>
parents: 222
diff changeset
37 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
38 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
39 end
345
7b5ef8b89268 Rungekutta bug fixes.
Jonatan Werpers <jonatan@werpers.com>
parents: 222
diff changeset
40 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
41 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
42 end
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
43
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
44 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
45 O = 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
46
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 obj.M = @(t)[
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 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
49 D(t), E(t);
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 ];
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 obj.C = @(t)[
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
52 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
53 S(t);
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 ];
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
55
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
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
60 obj.F = @(w,t)(obj.M(t)*w + obj.C(t));
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
61 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
62
e7e73173d44d time.Rungekutta4SecondOrder: Started adding code to allow timedependent D,E,S. Might be buggy.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
63 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
64 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
65 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
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 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
68 O = sparse(obj.m,obj.m);
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
69
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
70 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
71 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
72 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
73 ];
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.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
75 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
76 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 ];
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
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
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.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
84 end
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
85 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
86
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
87 function [v,t] = getV(obj)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
88 v = obj.w(1:end/2);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
89 t = obj.t;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
90 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
91
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
92 function [vt,t] = getVt(obj)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
93 vt = obj.w(end/2+1:end);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
94 t = obj.t;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
95 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
96
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
97 function obj = step(obj)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
98 obj.w = time.rk4.rungekutta_4(obj.w, obj.t, obj.k, obj.F);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
99 obj.t = obj.t + obj.k;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
100 obj.n = obj.n + 1;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
101 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
102 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
103
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
104
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
105 methods (Static)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
106 function k = getTimeStep(lambda)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
107 k = rk4.get_rk4_time_step(lambda);
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
108 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
109 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
110
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
111 end