Mercurial > repos > public > sbplib

comparison +time/+rk/ButcherTableau.m @ 990:1066bb31bc95 feature/timesteppers

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