Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_variable_6_2.m @ 321:5c9e5ba1c1ab feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clean up of d4_variable_6_2.m
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 26 Sep 2016 09:04:57 +0200
parents 99005a80b4c2
children c0cbffcf6513
comparison
equal deleted inserted replaced
320:3955ccffb7c8 321:5c9e5ba1c1ab
18 18
19 % H?r med 6 RP ist?llet f?r 8 f?r D4 operatorn, dock samma randderivator 19 % H?r med 6 RP ist?llet f?r 8 f?r D4 operatorn, dock samma randderivator
20 % Denna ?r noggrannare, och har 2a ordningens randdslutning och b?r ge 6te 20 % Denna ?r noggrannare, och har 2a ordningens randdslutning och b?r ge 6te
21 % ordningens konvergens. Hade dock ingen fri parameter att optimera 21 % ordningens konvergens. Hade dock ingen fri parameter att optimera
22 22
23 BP = 6;
24 if(m<2*BP)
25 error(['Operator requires at least ' num2str(2*BP) ' grid points']);
26 end
23 27
24 H = diag(ones(m,1),0); 28 % Norm
25 H(1:6,1:6) = [ 29 Hv = ones(m,1);
26 0.181e3/0.576e3 0 0 0 0 0; 30 Hv(1:6) = [0.181e3/0.576e3, 0.1343e4/0.960e3, 0.293e3/0.480e3, 0.1811e4/0.1440e4, 0.289e3/0.320e3, 0.65e2/0.64e2];
27 0 0.1343e4/0.960e3 0 0 0 0; 31 Hv(m-5:m) = rot90(Hv(1:6),2);
28 0 0 0.293e3/0.480e3 0 0 0; 32 Hv = h*Hv;
29 0 0 0 0.1811e4/0.1440e4 0 0; 33 H = spdiag(Hv, 0);
30 0 0 0 0 0.289e3/0.320e3 0; 34 HI = spdiag(1./Hv, 0);
31 0 0 0 0 0 0.65e2/0.64e2;
32 ];
33 H(m-5:m,m-5:m) = fliplr(flipud(H(1:6,1:6)));
34 35
35 36
36 e_1 = zeros(m,1);e_1(1) = 1; 37 % Boundary operators
37 e_m = zeros(m,1);e_m(m) = 1; 38 e_l = sparse(m,1);
39 e_l(1) = 1;
40 e_r = rot90(e_l, 2);
38 41
39 S_U = [-0.137e3/0.60e2 5 -5 0.10e2/0.3e1 -0.5e1/0.4e1 0.1e1/0.5e1;]/h; 42 d1_l = sparse(m,1);
40 S_1 = zeros(1,m); 43 d1_l(1:6) = [-0.137e3/0.60e2 5 -5 0.10e2/0.3e1 -0.5e1/0.4e1 0.1e1/0.5e1;]/h;
41 S_1(1:6) = S_U; 44 d1_r = -rot90(d1_l);
42 S_m = zeros(1,m);
43 S_m(m-5:m) = fliplr(-S_U);
44 45
45 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; 46 d2_l = sparse(m,1);
46 S2_1 = zeros(1,m); 47 d2_l(1:6) = [0.15e2/0.4e1 -0.77e2/0.6e1 0.107e3/0.6e1 -13 0.61e2/0.12e2 -0.5e1/0.6e1;]/h^2;
47 S2_1(1:6) = S2_U; 48 d2_r = rot90(d2_l, 2);
48 S2_m = zeros(1,m);
49 S2_m(m-5:m) = fliplr(S2_U);
50 49
51 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; 50 d3_l = sparse(m,1);
52 S3_1 = zeros(1,m); 51 d3_l(1:6) = [-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;
53 S3_1(1:6) = S3_U; 52 d3_r = -rot90(d3_l, 2);
54 S3_m = zeros(1,m);
55 S3_m(m-5:m) = fliplr(-S3_U);
56 53
57
58 H = h*H;
59 HI = inv(H);
60 54
61 % Fourth derivative, 1th order accurate at first 8 boundary points (still 55 % Fourth derivative, 1th order accurate at first 8 boundary points (still
62 % yield 5th order convergence if stable: for example u_tt = -u_xxxx 56 % yield 5th order convergence if stable: for example u_tt = -u_xxxx
63 57 stencil = [7/240, -2/5, 169/60, -122/15, 91/8, -122/15, 169/60, -2/5, 7/240];
64 m4 = 7/240; 58 diags = -4:4;
65 m3 = -2/5; 59 M4 = stripeMatrix(stencil, diags, m);
66 m2 = 169/60;
67 m1 = -122/15;
68 m0 = 91/8;
69
70 M4 = m4*(diag(ones(m-4,1),4)+diag(ones(m-4,1),-4))+m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0);
71
72 %M4 = (-1/6*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3) ) + 2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2)) -13/2*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1)) + 28/3*diag(ones(m,1),0));
73 60
74 M4_U = [ 61 M4_U = [
75 0.1009e4/0.192e3 -0.7657e4/0.480e3 0.9307e4/0.480e3 -0.509e3/0.40e2 0.4621e4/0.960e3 -0.25e2/0.32e2; 62 0.1009e4/0.192e3 -0.7657e4/0.480e3 0.9307e4/0.480e3 -0.509e3/0.40e2 0.4621e4/0.960e3 -0.25e2/0.32e2;
76 -0.7657e4/0.480e3 0.49513e5/0.960e3 -0.4007e4/0.60e2 0.21799e5/0.480e3 -0.8171e4/0.480e3 0.2657e4/0.960e3; 63 -0.7657e4/0.480e3 0.49513e5/0.960e3 -0.4007e4/0.60e2 0.21799e5/0.480e3 -0.8171e4/0.480e3 0.2657e4/0.960e3;
77 0.9307e4/0.480e3 -0.4007e4/0.60e2 0.1399e4/0.15e2 -0.2721e4/0.40e2 0.12703e5/0.480e3 -0.521e3/0.120e3; 64 0.9307e4/0.480e3 -0.4007e4/0.60e2 0.1399e4/0.15e2 -0.2721e4/0.40e2 0.12703e5/0.480e3 -0.521e3/0.120e3;
78 -0.509e3/0.40e2 0.21799e5/0.480e3 -0.2721e4/0.40e2 0.3349e4/0.60e2 -0.389e3/0.15e2 0.559e3/0.96e2; 65 -0.509e3/0.40e2 0.21799e5/0.480e3 -0.2721e4/0.40e2 0.3349e4/0.60e2 -0.389e3/0.15e2 0.559e3/0.96e2;
79 0.4621e4/0.960e3 -0.8171e4/0.480e3 0.12703e5/0.480e3 -0.389e3/0.15e2 0.17857e5/0.960e3 -0.1499e4/0.160e3; 66 0.4621e4/0.960e3 -0.8171e4/0.480e3 0.12703e5/0.480e3 -0.389e3/0.15e2 0.17857e5/0.960e3 -0.1499e4/0.160e3;
80 -0.25e2/0.32e2 0.2657e4/0.960e3 -0.521e3/0.120e3 0.559e3/0.96e2 -0.1499e4/0.160e3 0.2225e4/0.192e3; 67 -0.25e2/0.32e2 0.2657e4/0.960e3 -0.521e3/0.120e3 0.559e3/0.96e2 -0.1499e4/0.160e3 0.2225e4/0.192e3;
81 ]; 68 ];
82 69
70
83 M4(1:6,1:6) = M4_U; 71 M4(1:6,1:6) = M4_U;
72 M4(m-5:m,m-5:m) = rot90(M4_U, 2);
73 M4 = 1/h^3*M4;
84 74
85 M4(m-5:m,m-5:m) = flipud( fliplr( M4_U ) ); 75 D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r');
86 M4 = M4/h^3;
87
88 D4 = HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m);
89 end 76 end