Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_variable_6_min_boundary_points.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
20 20
21 H = diag(ones(m,1),0); 21 H = diag(ones(m,1),0);
22 H(1:6,1:6) = diag([13649/43200,12013/8640,2711/4320,5359/4320,7877/8640, 43801/43200]); 22 H(1:6,1:6) = diag([13649/43200,12013/8640,2711/4320,5359/4320,7877/8640, 43801/43200]);
23 H(m-5:m,m-5:m) = fliplr(flipud(diag([13649/43200,12013/8640, 2711/4320,5359/4320,7877/8640,43801/43200]))); 23 H(m-5:m,m-5:m) = fliplr(flipud(diag([13649/43200,12013/8640, 2711/4320,5359/4320,7877/8640,43801/43200])));
24 24
25 e_1 = zeros(m,1);e_1(1)=1;
26 e_m = zeros(m,1);e_m(m)=1;
27
28 S_U = [-25/12, 4, -3, 4/3, -1/4]/h;
29 S_1 = zeros(1,m);
30 S_1(1:5) = S_U;
31 S_m = zeros(1,m);
32 S_m(m-4:m) = fliplr(-S_U);
33
34 S2_U = [0.35e2/0.12e2 -0.26e2/0.3e1 0.19e2/0.2e1 -0.14e2/0.3e1 0.11e2/0.12e2;]/h^2;
35 S2_1 = zeros(1,m);
36 S2_1(1:5) = S2_U;
37 S2_m = zeros(1,m);
38 S2_m(m-4:m) = fliplr(S2_U);
39
40 S3_U = [-0.5e1/0.2e1 9 -12 7 -0.3e1/0.2e1;]/h^3;
41 S3_1 = zeros(1,m);
42 S3_1(1:5) = S3_U;
43 S3_m = zeros(1,m);
44 S3_m(m-4:m) = fliplr(-S3_U);
45
25 46
26 x1=0.70127127127127; 47 x1=0.70127127127127;
27 48
28 49
29 D1=(1/60*diag(ones(m-3,1),3)-9/60*diag(ones(m-2,1),2)+45/60*diag(ones(m-1,1),1)-45/60*diag(ones(m-1,1),-1)+9/60*diag(ones(m-2,1),-2)-1/60*diag(ones(m-3,1),-3)); 50 D1=(1/60*diag(ones(m-3,1),3)-9/60*diag(ones(m-2,1),2)+45/60*diag(ones(m-1,1),1)-45/60*diag(ones(m-1,1),-1)+9/60*diag(ones(m-2,1),-2)-1/60*diag(ones(m-3,1),-3));
36 34560/7877*x1-147127/47262, -129600/7877*x1+91715/7877, 172800/7877*x1-242111/15754, -86400/7877*x1+165041/23631, 0, 8640/7877*x1, -1296/7877, 144/7877, 0; 57 34560/7877*x1-147127/47262, -129600/7877*x1+91715/7877, 172800/7877*x1-242111/15754, -86400/7877*x1+165041/23631, 0, 8640/7877*x1, -1296/7877, 144/7877, 0;
37 -43200/43801*x1+89387/131403, 172800/43801*x1-240569/87602, -259200/43801*x1+182261/43801, 172800/43801*x1-710473/262806, -43200/43801*x1, 0, 32400/43801, -6480/43801, 720/43801 58 -43200/43801*x1+89387/131403, 172800/43801*x1-240569/87602, -259200/43801*x1+182261/43801, 172800/43801*x1-710473/262806, -43200/43801*x1, 0, 32400/43801, -6480/43801, 720/43801
38 ]; 59 ];
39 D1(m-5:m,m-8:m) = flipud( fliplr(-D1(1:6,1:9))); 60 D1(m-5:m,m-8:m) = flipud( fliplr(-D1(1:6,1:9)));
40 D1 = D1/h; 61 D1 = D1/h;
41
42 e_1 = zeros(m,1);e_1(1)=1;
43 e_m = zeros(m,1);e_m(m)=1;
44
45 S_U = [-25/12, 4, -3, 4/3, -1/4]/h;
46 S_1 = zeros(1,m);
47 S_1(1:5) = S_U;
48 S_m = zeros(1,m);
49 S_m(m-4:m) = fliplr(-S_U);
50
51
52 %DS=zeros(m,m);
53 %DS(1,1:5)=-[-25/12, 4, -3, 4/3, -1/4];
54 %DS(m,m-4:m)=fliplr(-[-25/12, 4, -3, 4/3, -1/4]);
55 %DS=diag(c)*DS/h;
56 62
57 63
58 H = h*H; 64 H = h*H;
59 HI = inv(H); 65 HI = inv(H);
60 66
95 101
96 M = M/h; 102 M = M/h;
97 103
98 D2 = HI*(-M-diag(c)*e_1*S_1+diag(c)*e_m*S_m); 104 D2 = HI*(-M-diag(c)*e_1*S_1+diag(c)*e_m*S_m);
99 105
100 S2_U = [0.35e2/0.12e2 -0.26e2/0.3e1 0.19e2/0.2e1 -0.14e2/0.3e1 0.11e2/0.12e2;]/h^2;
101 S2_1 = zeros(1,m);
102 S2_1(1:5) = S2_U;
103 S2_m = zeros(1,m);
104 S2_m(m-4:m) = fliplr(S2_U);
105
106
107
108
109 106
110 % Fourth derivative, 1th order accurate at first 8 boundary points (still 107 % Fourth derivative, 1th order accurate at first 8 boundary points (still
111 % yield 5th order convergence if stable: for example u_tt=-u_xxxx 108 % yield 5th order convergence if stable: for example u_tt=-u_xxxx
112 109
113 m4 = 7/240; 110 m4 = 7/240;
132 M4(1:6,1:6) = M4_U; 129 M4(1:6,1:6) = M4_U;
133 130
134 M4(m-5:m,m-5:m) = flipud( fliplr( M4_U ) ); 131 M4(m-5:m,m-5:m) = flipud( fliplr( M4_U ) );
135 M4 = M4/h^3; 132 M4 = M4/h^3;
136 133
137 S3_U = [-0.5e1/0.2e1 9 -12 7 -0.3e1/0.2e1;]/h^3;
138 S3_1 = zeros(1,m);
139 S3_1(1:5) = S3_U;
140 S3_m = zeros(1,m);
141 S3_m(m-4:m) = fliplr(-S3_U);
142
143 D4 = HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m); 134 D4 = HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m);
144 135
145 end 136 end