Mercurial > repos > public > sbplib

annotate +time/+rk/ButcherTableau.m @ 994:2f89959fb9f0 feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implement method to get gain from butcher tableu
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 09 Jan 2019 12:14:30 +0100
parents 1066bb31bc95
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
990
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
1 classdef ButcherTableau
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2 properties
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
3 a,b,c
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
4 end
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
5
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6 methods
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
7 % A ButcherTableau describes a specific rungekutta method where
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8 % y(n+1) = y(n) + dt*b(i)*k(i)
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9 % k(i) = F(t + c(i)*dt, Y(i))
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
10 % Y(i) = y(i) + dt*a(i,j)*k(j)
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
11 % where repeating indecies imply summation
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
12 function obj = ButcherTableau(a,b,c)
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
13 s = length(c);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
14 assertSize(a, [s,s]);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
15 assertLength(b, s);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
16
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
17 obj.a = a;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
18 obj.b = b;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
19 obj.c = c;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
20 end
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
21
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
22 function s = nStages(obj)
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
23 s = length(obj.c);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
24 end
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
25
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
26 function b = isExplicit(obj)
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
27 b = all(all(triu(obj.a)==0));
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
28 end
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
29
994
2f89959fb9f0 Implement method to get gain from butcher tableu
Jonatan Werpers <jonatan@werpers.com>
parents: 990
diff changeset
30 function g = testEquationGain(obj, z)
2f89959fb9f0 Implement method to get gain from butcher tableu
Jonatan Werpers <jonatan@werpers.com>
parents: 990
diff changeset
31 default_arg('z', sym('z'));
2f89959fb9f0 Implement method to get gain from butcher tableu
Jonatan Werpers <jonatan@werpers.com>
parents: 990
diff changeset
32 s = obj.nStages();
2f89959fb9f0 Implement method to get gain from butcher tableu
Jonatan Werpers <jonatan@werpers.com>
parents: 990
diff changeset
33
2f89959fb9f0 Implement method to get gain from butcher tableu
Jonatan Werpers <jonatan@werpers.com>
parents: 990
diff changeset
34 b = sym(obj.b);
2f89959fb9f0 Implement method to get gain from butcher tableu
Jonatan Werpers <jonatan@werpers.com>
parents: 990
diff changeset
35 A = sym(obj.a);
2f89959fb9f0 Implement method to get gain from butcher tableu
Jonatan Werpers <jonatan@werpers.com>
parents: 990
diff changeset
36 one = sym(ones(s,1));
2f89959fb9f0 Implement method to get gain from butcher tableu
Jonatan Werpers <jonatan@werpers.com>
parents: 990
diff changeset
37 I = sym(eye(s));
2f89959fb9f0 Implement method to get gain from butcher tableu
Jonatan Werpers <jonatan@werpers.com>
parents: 990
diff changeset
38
2f89959fb9f0 Implement method to get gain from butcher tableu
Jonatan Werpers <jonatan@werpers.com>
parents: 990
diff changeset
39 g = abs(1 + z*b*inv(I-z*A)*one);
2f89959fb9f0 Implement method to get gain from butcher tableu
Jonatan Werpers <jonatan@werpers.com>
parents: 990
diff changeset
40 end
2f89959fb9f0 Implement method to get gain from butcher tableu
Jonatan Werpers <jonatan@werpers.com>
parents: 990
diff changeset
41
990
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
42 % TBD: Add functions for checking accuracy, stability?
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
43 end
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
44
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
45 methods(Static)
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
46 % TVD (Total Variational Diminishing)
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
47 function bt = tvd_3()
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
48 a = [
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
49 0, 0, 0;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
50 1, 0, 0;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
51 1/4, 1/4, 0;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
52 ];
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
53 b = [1/6, 1/6, 2/3];
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
54 c = [0 1 1/2];
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
55
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
56 bt = time.rk.ButcherTableau(a,b,c);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
57 end
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
58
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
59 % Standard RK4
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
60 function bt = rk4()
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
61 a = [
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
62 0, 0, 0, 0;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
63 1/2, 0, 0, 0;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
64 0, 1/2, 0, 0;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
65 0, 0, 1, 0;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
66 ];
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
67
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
68 b = [1/6 1/3 1/3 1/6];
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
69 c = [0, 1/2, 1/2, 1];
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
70
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
71 bt = time.rk.ButcherTableau(a,b,c);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
72 end
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
73
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
74 % 3/8 RK4 (Kuttas method). Lower truncation error, more flops.
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
75 % Irreducible, unlike standard rk4.
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
76 function bt = rk4_3_8()
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
77 a = [
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
78 0, 0, 0, 0;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
79 1/3, 0, 0, 0;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
80 -1/3, 1, 0, 0;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
81 1, -1, 1, 0;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
82 ];
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
83
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
84 b = [1/8 3/8 3/8 1/8];
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
85 c = [0, 1/3, 2/3, 1];
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
86
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
87 bt = time.rk.ButcherTableau(a,b,c);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
88 end
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
89
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
90 % Runge-Kutta 6 from Alshina07
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
91 function bt = rk6()
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
92 a = zeros(7,7);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
93
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
94 a(2,1) = 4/7;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
95
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
96 a(3,1) = 115/112;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
97 a(3,2) = -5/16;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
98
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
99 a(4,1) = 589/630;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
100 a(4,2) = 5/18;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
101 a(4,3) = -16/45;
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
102
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
103 a(5,1) = 229/1200 - 29/6000*sqrt(5);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
104 a(5,2) = 119/240 - 187/1200*sqrt(5);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
105 a(5,3) = -14/75 + 34/375*sqrt(5);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
106 a(5,4) = -3/100*sqrt(5);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
107
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
108 a(6,1) = 71/2400 - 587/12000*sqrt(5);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
109 a(6,2) = 187/480 - 391/2400*sqrt(5);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
110 a(6,3) = -38/75 + 26/375*sqrt(5);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
111 a(6,4) = 27/80 - 3/400*sqrt(5);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
112 a(6,5) = (1+sqrt(5))/4
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
113
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
114 a(7,1) = -49/480 + 43/160*sqrt(5);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
115 a(7,2) = -425/96 + 51/32*sqrt(5);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
116 a(7,3) = 52/15 - 4/5*sqrt(5);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
117 a(7,4) = -27/16 + 3/16*sqrt(5);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
118 a(7,5) = 5/4 - 3/4*sqrt(5);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
119 a(7,6) = 5/2 - 1/2*sqrt(5);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
120
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
121 b = [1/12 0 0 0 5/12 5/12 1/12];
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
122 c = [0, 4/7, 5/7, 6/7, (5-sqrt(5))/10, (5+sqrt(5))/10, 1];
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
123
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
124 bt = time.rk.ButcherTableau(a,b,c);
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
125 end
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
126 end
1066bb31bc95 Create class for butcher tableau
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
127 end