Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_variable_6_3.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
46 S_1(1:6)=S_U; 46 S_1(1:6)=S_U;
47 S_m=zeros(1,m); 47 S_m=zeros(1,m);
48 S_m(m-5:m)=fliplr(-S_U); 48 S_m(m-5:m)=fliplr(-S_U);
49 49
50 50
51 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;
52 S2_1 = zeros(1,m);
53 S2_1(1:6) = S2_U;
54 S2_m = zeros(1,m);
55 S2_m(m-5:m) = fliplr(S2_U);
56
57
58 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;
59 S3_1 = zeros(1,m);
60 S3_1(1:6) = S3_U;
61 S3_m = zeros(1,m);
62 S3_m(m-5:m) = fliplr(-S3_U);
51 63
52 %DS=zeros(m,m); 64 %DS=zeros(m,m);
53 %DS(1,1:5)=-[-25/12, 4, -3, 4/3, -1/4]; 65 %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]); 66 %DS(m,m-4:m)=fliplr(-[-25/12, 4, -3, 4/3, -1/4]);
55 %DS=diag(c)*DS/h; 67 %DS=diag(c)*DS/h;
95 % 107 %
96 % M=M/h; 108 % M=M/h;
97 % 109 %
98 % D2=HI*(-M-diag(c)*e_1*S_1+diag(c)*e_m*S_m); 110 % D2=HI*(-M-diag(c)*e_1*S_1+diag(c)*e_m*S_m);
99 111
100 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;
101 S2_1 = zeros(1,m);
102 S2_1(1:6) = S2_U;
103 S2_m = zeros(1,m);
104 S2_m(m-5:m) = fliplr(S2_U);
105
106
107
108
109 112
110 % Fourth derivative, 1th order accurate at first 8 boundary points (still 113 % Fourth derivative, 1th order accurate at first 8 boundary points (still
111 % yield 5th order convergence if stable: for example u_tt=-u_xxxx 114 % yield 5th order convergence if stable: for example u_tt=-u_xxxx
112 115
113 m4 = 7/240; 116 m4 = 7/240;
133 M4(1:7,1:7) = M4_U; 136 M4(1:7,1:7) = M4_U;
134 137
135 M4(m-6:m,m-6:m) = flipud( fliplr( M4_U ) ); 138 M4(m-6:m,m-6:m) = flipud( fliplr( M4_U ) );
136 M4 = M4/h^3; 139 M4 = M4/h^3;
137 140
138 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;
139 S3_1 = zeros(1,m);
140 S3_1(1:6) = S3_U;
141 S3_m = zeros(1,m);
142 S3_m(m-5:m) = fliplr(-S3_U);
143
144 D4 = HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m); 141 D4 = HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m);
145 end 142 end