Mercurial > repos > public > sbplib

comparison +time/+rk/butcherTableau.m @ 916:d1c1615bd1a5 feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add comment in butcherTableau about rk4-3/8 being irreducible.
author Martin Almquist <malmquist@stanford.edu>
date Mon, 26 Nov 2018 16:08:54 -0800
parents 8732d6bd9890
children 0344fff87139
comparison
equal deleted inserted replaced
889:f5e14e5986b5 916:d1c1615bd1a5
17 c = [0 1 1/2]; 17 c = [0 1 1/2];
18 case "rk4" 18 case "rk4"
19 % Standard RK4 19 % Standard RK4
20 s = 4; 20 s = 4;
21 a = zeros(s,s-1); 21 a = zeros(s,s-1);
22 a(2,1) = 1/2; 22 a(2,1) = 1/2;
23 a(3,1) = 0; a(3,2) = 1/2; 23 a(3,1) = 0; a(3,2) = 1/2;
24 a(4,1) = 0; a(4,2) = 0; a(4,3) = 1; 24 a(4,1) = 0; a(4,2) = 0; a(4,3) = 1;
25 b = [1/6 1/3 1/3 1/6]; 25 b = [1/6 1/3 1/3 1/6];
26 c = [0, 1/2, 1/2, 1]; 26 c = [0, 1/2, 1/2, 1];
27 case "rk4-3/8" 27 case "rk4-3/8"
28 % 3/8 RK4 (Kuttas method). Lower truncation error, more flops 28 % 3/8 RK4 (Kuttas method). Lower truncation error, more flops.
29 % Irreducible, unlike standard rk4.
29 s = 4; 30 s = 4;
30 a = zeros(s,s-1); 31 a = zeros(s,s-1);
31 a(2,1) = 1/3; 32 a(2,1) = 1/3;
32 a(3,1) = -1/3; a(3,2) = 1; 33 a(3,1) = -1/3; a(3,2) = 1;
33 a(4,1) = 1; a(4,2) = -1; a(4,3) = 1; 34 a(4,1) = 1; a(4,2) = -1; a(4,3) = 1;
34 b = [1/8 3/8 3/8 1/8]; 35 b = [1/8 3/8 3/8 1/8];
35 c = [0, 1/3, 2/3, 1]; 36 c = [0, 1/3, 2/3, 1];
36 case "rk6" 37 case "rk6"
37 % Runge-Kutta 6 from Alshina07 38 % Runge-Kutta 6 from Alshina07
38 s = 7; 39 s = 7;
39 a = zeros(s,s-1); 40 a = zeros(s,s-1);
40 a(2,1) = 4/7; 41 a(2,1) = 4/7;
41 a(3,1) = 115/112; a(3,2) = -5/16; 42 a(3,1) = 115/112; a(3,2) = -5/16;
42 a(4,1) = 589/630; a(4,2) = 5/18; a(4,3) = -16/45; 43 a(4,1) = 589/630; a(4,2) = 5/18; a(4,3) = -16/45;
43 a(5,1) = 229/1200 - 29/6000*sqrt(5); a(5,2) = 119/240 - 187/1200*sqrt(5); a(5,3) = -14/75 + 34/375*sqrt(5); a(5,4) = -3/100*sqrt(5); 44 a(5,1) = 229/1200 - 29/6000*sqrt(5); a(5,2) = 119/240 - 187/1200*sqrt(5); a(5,3) = -14/75 + 34/375*sqrt(5); a(5,4) = -3/100*sqrt(5);
44 a(6,1) = 71/2400 - 587/12000*sqrt(5); a(6,2) = 187/480 - 391/2400*sqrt(5); a(6,3) = -38/75 + 26/375*sqrt(5); a(6,4) = 27/80 - 3/400*sqrt(5); a(6,5) = (1+sqrt(5))/4; 45 a(6,1) = 71/2400 - 587/12000*sqrt(5); a(6,2) = 187/480 - 391/2400*sqrt(5); a(6,3) = -38/75 + 26/375*sqrt(5); a(6,4) = 27/80 - 3/400*sqrt(5); a(6,5) = (1+sqrt(5))/4;
45 a(7,1) = -49/480 + 43/160*sqrt(5); a(7,2) = -425/96 + 51/32*sqrt(5); a(7,3) = 52/15 - 4/5*sqrt(5); a(7,4) = -27/16 + 3/16*sqrt(5); a(7,5) = 5/4 - 3/4*sqrt(5); a(7,6) = 5/2 - 1/2*sqrt(5); 46 a(7,1) = -49/480 + 43/160*sqrt(5); a(7,2) = -425/96 + 51/32*sqrt(5); a(7,3) = 52/15 - 4/5*sqrt(5); a(7,4) = -27/16 + 3/16*sqrt(5); a(7,5) = 5/4 - 3/4*sqrt(5); a(7,6) = 5/2 - 1/2*sqrt(5);
46 b = [1/12 0 0 0 5/12 5/12 1/12]; 47 b = [1/12 0 0 0 5/12 5/12 1/12];
47 c = [0, 4/7, 5/7, 6/7, (5-sqrt(5))/10, (5+sqrt(5))/10, 1]; 48 c = [0, 4/7, 5/7, 6/7, (5-sqrt(5))/10, (5+sqrt(5))/10, 1];
48 otherwise 49 otherwise
49 error('That Runge-Kutta method is not implemented', method) 50 error('That Runge-Kutta method is not implemented', method)
50 end 51 end