Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_variable_8_higher_boundary_order.m @ 316:203afa156f59 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Collected boundary operators.
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 23 Sep 2016 23:10:44 +0200
parents 9230c056a574
children 99005a80b4c2
comparison
equal deleted inserted replaced
315:297d2cbfbe15 316:203afa156f59
29 0 0 0 0 0 0 0 0.25641187e8/0.25401600e8; 29 0 0 0 0 0 0 0 0.25641187e8/0.25401600e8;
30 ]; 30 ];
31 31
32 H(m-7:m,m-7:m) = fliplr(flipud(H(1:8,1:8))); 32 H(m-7:m,m-7:m) = fliplr(flipud(H(1:8,1:8)));
33 33
34 e_1 = zeros(m,1);e_1(1) = 1; 34 e_1 = zeros(m,1);
35 e_m = zeros(m,1);e_m(m) = 1; 35 e_1(1) = 1;
36 e_m = zeros(m,1);
37 e_m(m) = 1;
36 38
37 S_U = [-0.49e2/0.20e2 6 -0.15e2/0.2e1 0.20e2/0.3e1 -0.15e2/0.4e1 0.6e1/0.5e1 -0.1e1/0.6e1]/h; 39 S_U = [-0.49e2/0.20e2 6 -0.15e2/0.2e1 0.20e2/0.3e1 -0.15e2/0.4e1 0.6e1/0.5e1 -0.1e1/0.6e1]/h;
38 S_1 = zeros(1,m); 40 S_1 = zeros(1,m);
39 S_1(1:7) = S_U; 41 S_1(1:7) = S_U;
40 S_m = zeros(1,m); 42 S_m = zeros(1,m);
41 S_m(m-6:m) = fliplr(-S_U); 43 S_m(m-6:m) = fliplr(-S_U);
44
45 S2_U = [0.203e3/0.45e2 -0.87e2/0.5e1 0.117e3/0.4e1 -0.254e3/0.9e1 0.33e2/0.2e1 -0.27e2/0.5e1 0.137e3/0.180e3]/h^2;
46 S2_1 = zeros(1,m);
47 S2_1(1:7) = S2_U;
48 S2_m = zeros(1,m);
49 S2_m(m-6:m) = fliplr(S2_U);
50
51 S3_U = [-0.49e2/0.8e1 29 -0.461e3/0.8e1 62 -0.307e3/0.8e1 13 -0.15e2/0.8e1]/h^3;
52 S3_1 = zeros(1,m);
53 S3_1(1:7) = S3_U;
54 S3_m = zeros(1,m);
55 S3_m(m-6:m) = fliplr(-S3_U);
42 56
43 H = h*H; 57 H = h*H;
44 HI = inv(H); 58 HI = inv(H);
45 59
46 % M = zeros(m,m); 60 % M = zeros(m,m);
79 % 93 %
80 % M = M/h; 94 % M = M/h;
81 % 95 %
82 % D2 = HI*(-M-diag(c)*e_1*S_1+diag(c)*e_m*S_m); 96 % D2 = HI*(-M-diag(c)*e_1*S_1+diag(c)*e_m*S_m);
83 97
84 S2_U = [0.203e3/0.45e2 -0.87e2/0.5e1 0.117e3/0.4e1 -0.254e3/0.9e1 0.33e2/0.2e1 -0.27e2/0.5e1 0.137e3/0.180e3]/h^2;
85 S2_1 = zeros(1,m);
86 S2_1(1:7) = S2_U;
87 S2_m = zeros(1,m);
88 S2_m(m-6:m) = fliplr(S2_U);
89
90
91 % Fourth derivative, 1th order accurate at first 8 boundary points (still 98 % Fourth derivative, 1th order accurate at first 8 boundary points (still
92 % yield 5th order convergence if stable: for example u_tt = -u_xxxx 99 % yield 5th order convergence if stable: for example u_tt = -u_xxxx
93 100
94 m5 = -0.41e2/0.7560e4; 101 m5 = -0.41e2/0.7560e4;
95 m4 = 0.1261e4/0.15120e5; 102 m4 = 0.1261e4/0.15120e5;
116 M4(1:8,1:8) = M4_U; 123 M4(1:8,1:8) = M4_U;
117 124
118 M4(m-7:m,m-7:m) = flipud( fliplr( M4_U ) ); 125 M4(m-7:m,m-7:m) = flipud( fliplr( M4_U ) );
119 M4 = M4/h^3; 126 M4 = M4/h^3;
120 127
121 S3_U = [-0.49e2/0.8e1 29 -0.461e3/0.8e1 62 -0.307e3/0.8e1 13 -0.15e2/0.8e1]/h^3;
122 S3_1 = zeros(1,m);
123 S3_1(1:7) = S3_U;
124 S3_m = zeros(1,m);
125 S3_m(m-6:m) = fliplr(-S3_U);
126
127 D4 = HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m); 128 D4 = HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m);
128 end 129 end