Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_variable_8_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
29 S_U = [-0.137e3/0.60e2 5 -5 0.10e2/0.3e1 -0.5e1/0.4e1 0.1e1/0.5e1;]/h; 29 S_U = [-0.137e3/0.60e2 5 -5 0.10e2/0.3e1 -0.5e1/0.4e1 0.1e1/0.5e1;]/h;
30 S_1 = zeros(1,m); 30 S_1 = zeros(1,m);
31 S_1(1:6) = S_U; 31 S_1(1:6) = S_U;
32 S_m = zeros(1,m); 32 S_m = zeros(1,m);
33 S_m(m-5:m) = fliplr(-S_U); 33 S_m(m-5:m) = fliplr(-S_U);
34
35 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;
36 S2_1 = zeros(1,m);
37 S2_1(1:6) = S2_U;
38 S2_m = zeros(1,m);
39 S2_m(m-5:m) = fliplr(S2_U);
40
41 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;
42 S3_1 = zeros(1,m);
43 S3_1(1:6) = S3_U;
44 S3_m = zeros(1,m);
45 S3_m(m-5:m) = fliplr(-S3_U);
46
47
34 48
35 H = h*H; 49 H = h*H;
36 HI = inv(H); 50 HI = inv(H);
37 51
38 52
72 % 86 %
73 % M=M/h; 87 % M=M/h;
74 % 88 %
75 % D2=HI*(-M-diag(c)*e_1*S_1+diag(c)*e_m*S_m); 89 % D2=HI*(-M-diag(c)*e_1*S_1+diag(c)*e_m*S_m);
76 90
77 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;
78 S2_1 = zeros(1,m);
79 S2_1(1:6) = S2_U;
80 S2_m = zeros(1,m);
81 S2_m(m-5:m) = fliplr(S2_U);
82
83 91
84 92
85 % Fourth derivative, 1th order accurate at first 8 boundary points (still 93 % Fourth derivative, 1th order accurate at first 8 boundary points (still
86 % yield 5th order convergence if stable: for example u_tt = -u_xxxx 94 % yield 5th order convergence if stable: for example u_tt = -u_xxxx
87 95
110 M4(1:8,1:8) = M4_U; 118 M4(1:8,1:8) = M4_U;
111 119
112 M4(m-7:m,m-7:m) = flipud( fliplr( M4_U ) ); 120 M4(m-7:m,m-7:m) = flipud( fliplr( M4_U ) );
113 M4 = M4/h^3; 121 M4 = M4/h^3;
114 122
115 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;
116 S3_1 = zeros(1,m);
117 S3_1(1:6) = S3_U;
118 S3_m = zeros(1,m);
119 S3_m(m-5:m) = fliplr(-S3_U);
120
121 D4=HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m); 123 D4=HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m);
122 end 124 end