Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_variable_6_2.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
40 S_1 = zeros(1,m); 40 S_1 = zeros(1,m);
41 S_1(1:6) = S_U; 41 S_1(1:6) = S_U;
42 S_m = zeros(1,m); 42 S_m = zeros(1,m);
43 S_m(m-5:m) = fliplr(-S_U); 43 S_m(m-5:m) = fliplr(-S_U);
44 44
45 S2_U = [0.15e2/0.4e1 -0.77e2/0.6e1 0.107e3/0.6e1 -13 0.61e2/0.12e2 -0.5e1/0.6e1;]/h^2;
46 S2_1 = zeros(1,m);
47 S2_1(1:6) = S2_U;
48 S2_m = zeros(1,m);
49 S2_m(m-5:m) = fliplr(S2_U);
45 50
46 51 S3_U = [-0.17e2/0.4e1 0.71e2/0.4e1 -0.59e2/0.2e1 0.49e2/0.2e1 -0.41e2/0.4e1 0.7e1/0.4e1;]/h^3;
47 %DS = zeros(m,m); 52 S3_1 = zeros(1,m);
48 %DS(1,1:5) = -[-25/12, 4, -3, 4/3, -1/4]; 53 S3_1(1:6) = S3_U;
49 %DS(m,m-4:m) = fliplr(-[-25/12, 4, -3, 4/3, -1/4]); 54 S3_m = zeros(1,m);
50 %DS = diag(c)*DS/h; 55 S3_m(m-5:m) = fliplr(-S3_U);
51 56
52 57
53 H = h*H; 58 H = h*H;
54 HI = inv(H); 59 HI = inv(H);
55 60
90 % 95 %
91 % M=M/h; 96 % M=M/h;
92 % 97 %
93 % D2=HI*(-M-diag(c)*e_1*S_1+diag(c)*e_m*S_m); 98 % D2=HI*(-M-diag(c)*e_1*S_1+diag(c)*e_m*S_m);
94 99
95 S2_U = [0.15e2/0.4e1 -0.77e2/0.6e1 0.107e3/0.6e1 -13 0.61e2/0.12e2 -0.5e1/0.6e1;]/h^2;
96 S2_1 = zeros(1,m);
97 S2_1(1:6) = S2_U;
98 S2_m = zeros(1,m);
99 S2_m(m-5:m) = fliplr(S2_U);
100
101
102
103
104 100
105 % Fourth derivative, 1th order accurate at first 8 boundary points (still 101 % Fourth derivative, 1th order accurate at first 8 boundary points (still
106 % yield 5th order convergence if stable: for example u_tt = -u_xxxx 102 % yield 5th order convergence if stable: for example u_tt = -u_xxxx
107 103
108 m4 = 7/240; 104 m4 = 7/240;
127 M4(1:6,1:6) = M4_U; 123 M4(1:6,1:6) = M4_U;
128 124
129 M4(m-5:m,m-5:m) = flipud( fliplr( M4_U ) ); 125 M4(m-5:m,m-5:m) = flipud( fliplr( M4_U ) );
130 M4 = M4/h^3; 126 M4 = M4/h^3;
131 127
132 S3_U = [-0.17e2/0.4e1 0.71e2/0.4e1 -0.59e2/0.2e1 0.49e2/0.2e1 -0.41e2/0.4e1 0.7e1/0.4e1;]/h^3;
133 S3_1 = zeros(1,m);
134 S3_1(1:6) = S3_U;
135 S3_m = zeros(1,m);
136 S3_m(m-5:m) = fliplr(-S3_U);
137
138 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);
139
140
141 end 129 end