Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d2_variable_periodic_12.m @ 1337:bf2554f1825d feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add periodic D1 and D2 operators for orders 8,10,12
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 13 May 2022 13:28:10 +0200
parents
children
comparison
equal deleted inserted replaced
1336:0666629aa183 1337:bf2554f1825d
1 function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_12(m,h)
2 % m = number of unique grid points, i.e. h = L/m;
3
4 if(m<13)
5 error(['Operator requires at least ' num2str(8) ' grid points']);
6 end
7
8 % Norm
9 Hv = ones(m,1);
10 Hv = h*Hv;
11 H = spdiag(Hv, 0);
12 HI = spdiag(1./Hv, 0);
13
14
15 % Dummy boundary operators
16 e_l = sparse(m,1);
17 e_r = rot90(e_l, 2);
18
19 d1_l = sparse(m,1);
20 d1_r = -rot90(d1_l, 2);
21
22
23 % D1 operator
24 diags = -6:6;
25 stencil = [1/5544, -1/385, 1/56, -5/63, 15/56, -6/7, 0, 6/7, -15/56, 5/63, -1/56, 1/385, -1/5544];
26 D1 = stripeMatrixPeriodic(stencil, diags, m);
27 D1 = D1/h;
28
29 % Undivided differences
30 diags = -4:3;
31 stencil = [-1 7 -21 35 -35 21 -7 1];
32 DD_7 = stripeMatrixPeriodic(stencil, diags, m);
33
34 diags = -4:4;
35 stencil = [1 -8 28 -56 70 -56 28 -8 1];
36 DD_8 = stripeMatrixPeriodic(stencil, diags, m);
37
38 diags = -5:4;
39 stencil = [-1 9 -36 84 -126 126 -84 36 -9 1];
40 DD_9 = stripeMatrixPeriodic(stencil, diags, m);
41
42 diags = -5:5;
43 stencil = [1 -10 45 -120 210 -252 210 -120 45 -10 1];
44 DD_10 = stripeMatrixPeriodic(stencil, diags, m);
45
46 diags = -6:5;
47 stencil = [-1, 11, -55, 165, -330, 462, -462, 330, -165, 55, -11, 1];
48 DD_11 = stripeMatrixPeriodic(stencil, diags, m);
49
50 diags = -6:6;
51 stencil = [1, -12, 66, -220, 495, -792, 924, -792, 495, -220, 66,-12, 1];
52 DD_12 = stripeMatrixPeriodic(stencil, diags, m);
53
54 % D2 operator
55 function D2 = D2_fun(c)
56
57 diags = -1:0;
58 stencil = [1/2, 1/2];
59 C2 = stripeMatrixPeriodic(stencil, diags, m);
60
61 diags = -1:1;
62 stencil = [5/18, 4/9, 5/18];
63 C3 = stripeMatrixPeriodic(stencil, diags, m);
64
65 diags = -2:1;
66 stencil = [5/28, 9/28, 9/28, 5/28];
67 C4 = stripeMatrixPeriodic(stencil, diags, m);
68
69 diags = -2:2;
70 stencil = [1/7, 8/35, 9/35, 8/35, 1/7];
71 C5 = stripeMatrixPeriodic(stencil, diags, m);
72
73 diags = -3:2;
74 stencil = [1/6, 1/6, 1/6, 1/6, 1/6, 1/6];
75 C6 = stripeMatrixPeriodic(stencil, diags, m);
76
77 C1 = sparse(diag(c));
78 C2 = sparse(diag(C2 * c));
79 C3 = sparse(diag(C3 * c));
80 C4 = sparse(diag(C4 * c));
81
82 % Remainder term added to wide second derivative operator
83 R = (1/30735936 / h) * transpose(DD_12) * C1 * DD_12 + (1/6403320 / h) * transpose(DD_11) * C2 * DD_11 + (1/1293600 / h) * transpose(DD_10) * C3 * DD_10 + (1/249480 / h) * transpose(DD_9) * C4 * DD_9 + (1/44352 / h) * transpose(DD_8) * C5 * DD_8 + (1/6468 / h) * transpose(DD_7) * C6 * DD_7;
84 D2 = D1 * C1 * D1 - H \ R;
85 end
86 D2 = @D2_fun;
87
88
89 end