Mercurial > repos > public > sbplib

annotate +sbp/+implementations/d2_variable_periodic_10.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
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1337
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
1 function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_10(m,h)
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
2 % m = number of unique grid points, i.e. h = L/m;
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
3
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
4 if(m<11)
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
5 error(['Operator requires at least ' num2str(8) ' grid points']);
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
6 end
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
7
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
8 % Norm
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
9 Hv = ones(m,1);
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
10 Hv = h*Hv;
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
11 H = spdiag(Hv, 0);
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
12 HI = spdiag(1./Hv, 0);
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
13
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
14
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
15 % Dummy boundary operators
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
16 e_l = sparse(m,1);
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
17 e_r = rot90(e_l, 2);
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
18
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
19 d1_l = sparse(m,1);
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
20 d1_r = -rot90(d1_l, 2);
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
21
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
22
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
23 % D1 operator
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
24 diags = -5:5;
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
25 stencil = [-1/1260, 5/504, -5/84, 5/21, -5/6, 0, 5/6, -5/21, 5/84, -5/504, 1/1260];
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
26 D1 = stripeMatrixPeriodic(stencil, diags, m);
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
27 D1 = D1/h;
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
28
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
29 % Undivided differences
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
30 diags = -3:3;
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
31 stencil = [1 -6 15 -20 15 -6 1];
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
32 DD_6 = stripeMatrixPeriodic(stencil, diags, m);
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
33
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
34 diags = -4:3;
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
35 stencil = [-1 7 -21 35 -35 21 -7 1];
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
36 DD_7 = stripeMatrixPeriodic(stencil, diags, m);
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
37
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
38 diags = -4:4;
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
39 stencil = [1 -8 28 -56 70 -56 28 -8 1];
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
40 DD_8 = stripeMatrixPeriodic(stencil, diags, m);
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
41
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
42 diags = -5:4;
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
43 stencil = [-1 9 -36 84 -126 126 -84 36 -9 1];
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
44 DD_9 = stripeMatrixPeriodic(stencil, diags, m);
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
45
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
46 diags = -5:5;
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
47 stencil = [1 -10 45 -120 210 -252 210 -120 45 -10 1];
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
48 DD_10 = stripeMatrixPeriodic(stencil, diags, m);
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
49 % D2 operator
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
50 function D2 = D2_fun(c)
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
51
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
52 diags = -1:0;
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
53 stencil = [1/2, 1/2];
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
54 C2 = stripeMatrixPeriodic(stencil, diags, m);
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
55
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
56 diags = -1:1;
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
57 stencil = [2/7, 3/7, 2/7];
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
58 C3 = stripeMatrixPeriodic(stencil, diags, m);
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
59
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
60 diags = -2:1;
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
61 stencil = [1/5, 3/10, 1/5, 3/10];
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
62 C4 = stripeMatrixPeriodic(stencil, diags, m);
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
63
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
64 diags = -2:2;
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
65 stencil = [1/5, 1/5, 1/5, 1/5, 1/5];
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
66 C5 = stripeMatrixPeriodic(stencil, diags, m);
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
67
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
68 C1 = sparse(diag(c));
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
69 C2 = sparse(diag(C2 * c));
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
70 C3 = sparse(diag(C3 * c));
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
71 C4 = sparse(diag(C4 * c));
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
72
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
73 % Remainder term added to wide second derivative operator
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
74 R = (1/1587600 / h) * transpose(DD_10) * C1 * DD_10 + (1/317520 / h) * transpose(DD_9) * C2 * DD_9 + (1/60480 / h) * transpose(DD_8) * C3 * DD_8 + (1/10584 / h) * transpose(DD_7) * C4 * DD_7 + (1/1512 / h) * transpose(DD_6) * C5 * DD_6;
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
75 D2 = D1 * C1 * D1 - H \ R;
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
76 end
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
77 D2 = @D2_fun;
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
78
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
79
bf2554f1825d Add periodic D1 and D2 operators for orders 8,10,12
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
diff changeset
80 end