Mercurial > repos > public > sbplib
comparison +sbp/+implementations/d4_variable_4.m @ 314:88584b0cfba1 feature/beams
Corrections and clean up order 4
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Fri, 23 Sep 2016 22:55:30 +0200 |
parents | 9230c056a574 |
children | 297d2cbfbe15 |
comparison
equal
deleted
inserted
replaced
313:52b4cdf27633 | 314:88584b0cfba1 |
---|---|
11 BP = 6; | 11 BP = 6; |
12 if(m<2*BP) | 12 if(m<2*BP) |
13 error(['Operator requires at least ' num2str(2*BP) ' grid points']); | 13 error(['Operator requires at least ' num2str(2*BP) ' grid points']); |
14 end | 14 end |
15 | 15 |
16 H = speye(m,m); | 16 % Norm |
17 H(1:4,1:4) = diag([17/48 59/48 43/48 49/48]); | 17 Hv = ones(m,1); |
18 H(m-3:m,m-3:m) = rot90(diag([17/48 59/48 43/48 49/48]),2); | 18 Hv(1:4) = [17/48 59/48 43/48 49/48]; |
19 H = H*h; | 19 Hv(m-3:m) = rot90(Hv(1:4),2); |
20 HI = inv(H); | 20 Hv = h*Hv; |
21 HI = sparse(HI); | 21 H = spdiag(Hv, 0); |
22 HI = spdiag(1./Hv, 0); | |
22 | 23 |
23 | 24 |
24 % Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2)); | 25 % Boundary operators |
25 e = ones(m,1); | 26 e_l = sparse(m,1); |
26 % Q=spdiags([e -8*e 0*e 8*e -e], -2:2, m, m)/12; | 27 e_l(1) = 1; |
27 % Q_U = [0 0.59e2/0.96e2 -0.1e1/0.12e2 -0.1e1/0.32e2; -0.59e2/0.96e2 0 0.59e2/0.96e2 0; 0.1e1/0.12e2 -0.59e2/0.96e2 0 0.59e2/0.96e2; 0.1e1/0.32e2 0 -0.59e2/0.96e2 0;]; | 28 e_r = rot90(e_l, 2); |
28 % Q(1:4,1:4)=Q_U; | |
29 % Q(m-3:m,m-3:m)=rot90( -Q_U(1:4,1:4) ,2 ); | |
30 | 29 |
31 e_1 = sparse(m,1); | 30 d1_l = sparse(m,1); |
32 e_1(1) = 1; | 31 d1_l(1:4) = 1/h*[-11/6 3 -3/2 1/3]; |
33 e_m = sparse(m,1); | 32 d1_r = -rot90(d1_l); |
34 e_m(m) = 1; | |
35 | 33 |
36 % D1=HI*(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; | 34 d2_l = sparse(m,1); |
35 d2_l(1:4) = 1/h^2*[2 -5 4 -1]; | |
36 d2_r = rot90(d2_l, 2); | |
37 | 37 |
38 M_U = [ | 38 d3_l = sparse(m,1); |
39 0.9e1/0.8e1 -0.59e2/0.48e2 0.1e1/0.12e2 0.1e1/0.48e2; | 39 d3_l(1:4) = 1/h^3*[-1 3 -3 1]; |
40 -0.59e2/0.48e2 0.59e2/0.24e2 -0.59e2/0.48e2 0; | 40 d3_r = -rot90(d3_l, 2); |
41 0.1e1/0.12e2 -0.59e2/0.48e2 0.55e2/0.24e2 -0.59e2/0.48e2; | |
42 0.1e1/0.48e2 0 -0.59e2/0.48e2 0.59e2/0.24e2; | |
43 ]; | |
44 % M=-(-1/12*diag(ones(m-2,1),2)+16/12*diag(ones(m-1,1),1)+16/12*diag(ones(m-1,1),-1)-1/12*diag(ones(m-2,1),-2)-30/12*diag(ones(m,1),0)); | |
45 M = -spdiags([-e 16*e -30*e 16*e -e], -2:2, m, m)/12; | |
46 | |
47 M(1:4,1:4) = M_U; | |
48 | |
49 M(m-3:m,m-3:m) = rot90( M_U ,2 ); | |
50 M=M/h; | |
51 | |
52 S_U=[-0.11e2/0.6e1 3 -0.3e1/0.2e1 0.1e1/0.3e1;]/h; | |
53 S_1=sparse(1,m); | |
54 S_1(1:4)=S_U; | |
55 S_m=sparse(1,m); | |
56 S_m(m-3:m)=fliplr(-S_U); | |
57 S_1 = S_1'; | |
58 S_m = S_m'; | |
59 | 41 |
60 | 42 |
43 % First derivative SBP operator, | |
44 stencil = [1/12 -2/3 0 2/3 -1/12]; | |
45 diags = [-1 0 1]; | |
46 | |
47 Q_U = [ | |
48 0 0.59e2/0.96e2 -0.1e1/0.12e2 -0.1e1/0.32e2; | |
49 -0.59e2/0.96e2 0 0.59e2/0.96e2 0; | |
50 0.1e1/0.12e2 -0.59e2/0.96e2 0 0.59e2/0.96e2; | |
51 0.1e1/0.32e2 0 -0.59e2/0.96e2 0; | |
52 ]; | |
53 | |
54 Q = stripeMatrix(stencil, diags, m); | |
55 Q(1:4,1:4)=Q_U; | |
56 Q(m-3:m,m-3:m) = -rot90(Q_U, 2); | |
57 | |
58 D1 = HI*(Q - 1/2*e_l*e_l' + 1/2*e_r*e_r'); | |
59 | |
60 | |
61 % Second derivative | |
61 M=sparse(m,m); | 62 M=sparse(m,m); |
62 e_1 = sparse(e_1); | |
63 e_m = sparse(e_m); | |
64 S_1 = sparse(S_1); | |
65 S_m = sparse(S_m); | |
66 | 63 |
67 scheme_width = 5; | 64 scheme_width = 5; |
68 scheme_radius = (scheme_width-1)/2; | 65 scheme_radius = (scheme_width-1)/2; |
69 r = (1+scheme_radius):(m-scheme_radius); | 66 r = (1+scheme_radius):(m-scheme_radius); |
70 | 67 |
107 % M(M_diag_ind) = [Mm2 Mm1 M0 Mp1 Mp2]; % This is slightly faster | 104 % M(M_diag_ind) = [Mm2 Mm1 M0 Mp1 Mp2]; % This is slightly faster |
108 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% | 105 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
109 % Kan man skriva det som en multiplikation av en 3-dim matris? | 106 % Kan man skriva det som en multiplikation av en 3-dim matris? |
110 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% | 107 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
111 | 108 |
112 | |
113 | |
114 | |
115 M(1:6,1:6) = [ | 109 M(1:6,1:6) = [ |
116 0.12e2/0.17e2 * c(1) + 0.59e2/0.192e3 * c(2) + 0.27010400129e11/0.345067064608e12 * c(3) + 0.69462376031e11/0.2070402387648e13 * c(4) -0.59e2/0.68e2 * c(1) - 0.6025413881e10/0.21126554976e11 * c(3) - 0.537416663e9/0.7042184992e10 * c(4) 0.2e1/0.17e2 * c(1) - 0.59e2/0.192e3 * c(2) + 0.213318005e9/0.16049630912e11 * c(4) + 0.2083938599e10/0.8024815456e10 * c(3) 0.3e1/0.68e2 * c(1) - 0.1244724001e10/0.21126554976e11 * c(3) + 0.752806667e9/0.21126554976e11 * c(4) 0.49579087e8/0.10149031312e11 * c(3) - 0.49579087e8/0.10149031312e11 * c(4) -c(4)/0.784e3 + c(3)/0.784e3; | 110 0.12e2/0.17e2 * c(1) + 0.59e2/0.192e3 * c(2) + 0.27010400129e11/0.345067064608e12 * c(3) + 0.69462376031e11/0.2070402387648e13 * c(4) -0.59e2/0.68e2 * c(1) - 0.6025413881e10/0.21126554976e11 * c(3) - 0.537416663e9/0.7042184992e10 * c(4) 0.2e1/0.17e2 * c(1) - 0.59e2/0.192e3 * c(2) + 0.213318005e9/0.16049630912e11 * c(4) + 0.2083938599e10/0.8024815456e10 * c(3) 0.3e1/0.68e2 * c(1) - 0.1244724001e10/0.21126554976e11 * c(3) + 0.752806667e9/0.21126554976e11 * c(4) 0.49579087e8/0.10149031312e11 * c(3) - 0.49579087e8/0.10149031312e11 * c(4) -c(4)/0.784e3 + c(3)/0.784e3; |
117 -0.59e2/0.68e2 * c(1) - 0.6025413881e10/0.21126554976e11 * c(3) - 0.537416663e9/0.7042184992e10 * c(4) 0.3481e4/0.3264e4 * c(1) + 0.9258282831623875e16/0.7669235228057664e16 * c(3) + 0.236024329996203e15/0.1278205871342944e16 * c(4) -0.59e2/0.408e3 * c(1) - 0.29294615794607e14/0.29725717938208e14 * c(3) - 0.2944673881023e13/0.29725717938208e14 * c(4) -0.59e2/0.1088e4 * c(1) + 0.260297319232891e15/0.2556411742685888e16 * c(3) - 0.60834186813841e14/0.1278205871342944e16 * c(4) -0.1328188692663e13/0.37594290333616e14 * c(3) + 0.1328188692663e13/0.37594290333616e14 * c(4) -0.8673e4/0.2904112e7 * c(3) + 0.8673e4/0.2904112e7 * c(4); | 111 -0.59e2/0.68e2 * c(1) - 0.6025413881e10/0.21126554976e11 * c(3) - 0.537416663e9/0.7042184992e10 * c(4) 0.3481e4/0.3264e4 * c(1) + 0.9258282831623875e16/0.7669235228057664e16 * c(3) + 0.236024329996203e15/0.1278205871342944e16 * c(4) -0.59e2/0.408e3 * c(1) - 0.29294615794607e14/0.29725717938208e14 * c(3) - 0.2944673881023e13/0.29725717938208e14 * c(4) -0.59e2/0.1088e4 * c(1) + 0.260297319232891e15/0.2556411742685888e16 * c(3) - 0.60834186813841e14/0.1278205871342944e16 * c(4) -0.1328188692663e13/0.37594290333616e14 * c(3) + 0.1328188692663e13/0.37594290333616e14 * c(4) -0.8673e4/0.2904112e7 * c(3) + 0.8673e4/0.2904112e7 * c(4); |
118 0.2e1/0.17e2 * c(1) - 0.59e2/0.192e3 * c(2) + 0.213318005e9/0.16049630912e11 * c(4) + 0.2083938599e10/0.8024815456e10 * c(3) -0.59e2/0.408e3 * c(1) - 0.29294615794607e14/0.29725717938208e14 * c(3) - 0.2944673881023e13/0.29725717938208e14 * c(4) c(1)/0.51e2 + 0.59e2/0.192e3 * c(2) + 0.13777050223300597e17/0.26218083221499456e17 * c(4) + 0.564461e6/0.13384296e8 * c(5) + 0.378288882302546512209e21/0.270764341349677687456e21 * c(3) c(1)/0.136e3 - 0.125059e6/0.743572e6 * c(5) - 0.4836340090442187227e19/0.5525802884687299744e19 * c(3) - 0.17220493277981e14/0.89177153814624e14 * c(4) -0.10532412077335e14/0.42840005263888e14 * c(4) + 0.1613976761032884305e19/0.7963657098519931984e19 * c(3) + 0.564461e6/0.4461432e7 * c(5) -0.960119e6/0.1280713392e10 * c(4) - 0.3391e4/0.6692148e7 * c(5) + 0.33235054191e11/0.26452850508784e14 * c(3); | 112 0.2e1/0.17e2 * c(1) - 0.59e2/0.192e3 * c(2) + 0.213318005e9/0.16049630912e11 * c(4) + 0.2083938599e10/0.8024815456e10 * c(3) -0.59e2/0.408e3 * c(1) - 0.29294615794607e14/0.29725717938208e14 * c(3) - 0.2944673881023e13/0.29725717938208e14 * c(4) c(1)/0.51e2 + 0.59e2/0.192e3 * c(2) + 0.13777050223300597e17/0.26218083221499456e17 * c(4) + 0.564461e6/0.13384296e8 * c(5) + 0.378288882302546512209e21/0.270764341349677687456e21 * c(3) c(1)/0.136e3 - 0.125059e6/0.743572e6 * c(5) - 0.4836340090442187227e19/0.5525802884687299744e19 * c(3) - 0.17220493277981e14/0.89177153814624e14 * c(4) -0.10532412077335e14/0.42840005263888e14 * c(4) + 0.1613976761032884305e19/0.7963657098519931984e19 * c(3) + 0.564461e6/0.4461432e7 * c(5) -0.960119e6/0.1280713392e10 * c(4) - 0.3391e4/0.6692148e7 * c(5) + 0.33235054191e11/0.26452850508784e14 * c(3); |
119 0.3e1/0.68e2 * c(1) - 0.1244724001e10/0.21126554976e11 * c(3) + 0.752806667e9/0.21126554976e11 * c(4) -0.59e2/0.1088e4 * c(1) + 0.260297319232891e15/0.2556411742685888e16 * c(3) - 0.60834186813841e14/0.1278205871342944e16 * c(4) c(1)/0.136e3 - 0.125059e6/0.743572e6 * c(5) - 0.4836340090442187227e19/0.5525802884687299744e19 * c(3) - 0.17220493277981e14/0.89177153814624e14 * c(4) 0.3e1/0.1088e4 * c(1) + 0.507284006600757858213e21/0.475219048083107777984e21 * c(3) + 0.1869103e7/0.2230716e7 * c(5) + c(6)/0.24e2 + 0.1950062198436997e16/0.3834617614028832e16 * c(4) -0.4959271814984644613e19/0.20965546238960637264e20 * c(3) - c(6)/0.6e1 - 0.15998714909649e14/0.37594290333616e14 * c(4) - 0.375177e6/0.743572e6 * c(5) -0.368395e6/0.2230716e7 * c(5) + 0.752806667e9/0.539854092016e12 * c(3) + 0.1063649e7/0.8712336e7 * c(4) + c(6)/0.8e1; | 113 0.3e1/0.68e2 * c(1) - 0.1244724001e10/0.21126554976e11 * c(3) + 0.752806667e9/0.21126554976e11 * c(4) -0.59e2/0.1088e4 * c(1) + 0.260297319232891e15/0.2556411742685888e16 * c(3) - 0.60834186813841e14/0.1278205871342944e16 * c(4) c(1)/0.136e3 - 0.125059e6/0.743572e6 * c(5) - 0.4836340090442187227e19/0.5525802884687299744e19 * c(3) - 0.17220493277981e14/0.89177153814624e14 * c(4) 0.3e1/0.1088e4 * c(1) + 0.507284006600757858213e21/0.475219048083107777984e21 * c(3) + 0.1869103e7/0.2230716e7 * c(5) + c(6)/0.24e2 + 0.1950062198436997e16/0.3834617614028832e16 * c(4) -0.4959271814984644613e19/0.20965546238960637264e20 * c(3) - c(6)/0.6e1 - 0.15998714909649e14/0.37594290333616e14 * c(4) - 0.375177e6/0.743572e6 * c(5) -0.368395e6/0.2230716e7 * c(5) + 0.752806667e9/0.539854092016e12 * c(3) + 0.1063649e7/0.8712336e7 * c(4) + c(6)/0.8e1; |
134 M = M/h; | 128 M = M/h; |
135 D2 = HI*(-M-c(1)*e_1*S_1'+c(m)*e_m*S_m'); | 129 D2 = HI*(-M-c(1)*e_1*S_1'+c(m)*e_m*S_m'); |
136 end | 130 end |
137 D2 = @D2_fun; | 131 D2 = @D2_fun; |
138 | 132 |
139 S2_U=[2 -5 4 -1;]/h^2; | |
140 S2_1=sparse(1,m); | |
141 S2_1(1:4)=S2_U; | |
142 S2_m=sparse(1,m); | |
143 S2_m(m-3:m)=fliplr(S2_U); | |
144 S2_1 = S2_1'; | |
145 S2_m = S2_m'; | |
146 | 133 |
147 m3 = -1/6; | 134 % Fourth derivative |
148 m2 = 2; | 135 stencil = [-1/6,2,-13/2, 28/3,-13/2,2,-1/6]; |
149 m1 = -13/2; | 136 diags = -3:3; |
150 m0 = 28/3; | |
151 % M4=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); | |
152 stencil = [m3,m2,m1,m0,m1,m2,m3]; | |
153 d = (length(stencil)-1)/2; | |
154 diags = -d:d; | |
155 M4 = stripeMatrix(stencil, diags, m); | 137 M4 = stripeMatrix(stencil, diags, m); |
156 | |
157 %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)); | |
158 | 138 |
159 M4_U = [ | 139 M4_U = [ |
160 0.5762947e7/0.2316384e7 -0.6374287e7/0.1158192e7 0.573947e6/0.165456e6 -0.124637e6/0.289548e6 0.67979e5/0.2316384e7 -0.60257e5/0.1158192e7; | 140 0.5762947e7/0.2316384e7 -0.6374287e7/0.1158192e7 0.573947e6/0.165456e6 -0.124637e6/0.289548e6 0.67979e5/0.2316384e7 -0.60257e5/0.1158192e7; |
161 -0.6374287e7/0.1158192e7 0.30392389e8/0.2316384e7 -0.2735053e7/0.289548e6 0.273109e6/0.165456e6 0.83767e5/0.1158192e7 0.245549e6/0.2316384e7; | 141 -0.6374287e7/0.1158192e7 0.30392389e8/0.2316384e7 -0.2735053e7/0.289548e6 0.273109e6/0.165456e6 0.83767e5/0.1158192e7 0.245549e6/0.2316384e7; |
162 0.573947e6/0.165456e6 -0.2735053e7/0.289548e6 0.5266855e7/0.579096e6 -0.1099715e7/0.289548e6 0.869293e6/0.1158192e7 -0.10195e5/0.144774e6; | 142 0.573947e6/0.165456e6 -0.2735053e7/0.289548e6 0.5266855e7/0.579096e6 -0.1099715e7/0.289548e6 0.869293e6/0.1158192e7 -0.10195e5/0.144774e6; |
164 0.67979e5/0.2316384e7 0.83767e5/0.1158192e7 0.869293e6/0.1158192e7 -0.324229e6/0.72387e5 0.2626501e7/0.330912e6 -0.7115491e7/0.1158192e7; | 144 0.67979e5/0.2316384e7 0.83767e5/0.1158192e7 0.869293e6/0.1158192e7 -0.324229e6/0.72387e5 0.2626501e7/0.330912e6 -0.7115491e7/0.1158192e7; |
165 -0.60257e5/0.1158192e7 0.245549e6/0.2316384e7 -0.10195e5/0.144774e6 0.1847891e7/0.1158192e7 -0.7115491e7/0.1158192e7 0.21383077e8/0.2316384e7; | 145 -0.60257e5/0.1158192e7 0.245549e6/0.2316384e7 -0.10195e5/0.144774e6 0.1847891e7/0.1158192e7 -0.7115491e7/0.1158192e7 0.21383077e8/0.2316384e7; |
166 ]; | 146 ]; |
167 | 147 |
168 M4(1:6,1:6) = M4_U; | 148 M4(1:6,1:6) = M4_U; |
149 M4(m-5:m,m-5:m) = rot90(M4_U, 2); | |
150 M4 = 1/h^3*M4; | |
169 | 151 |
170 M4(m-5:m,m-5:m) = rot90( M4_U ,2 ); | 152 D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r'); |
171 M4 = M4/h^3; | |
172 | |
173 S3_U = [-1 3 -3 1;]/h^3; | |
174 S3_1 = sparse(1,m); | |
175 S3_1(1:4)=S3_U; | |
176 S3_m = sparse(1,m); | |
177 S3_m(m-3:m) = fliplr(-S3_U); | |
178 S3_1 = S3_1'; | |
179 S3_m = S3_m'; | |
180 | |
181 D4=HI*(M4-e_1*S3_1'+e_m*S3_m' + S_1*S2_1'-S_m*S2_m'); | |
182 | |
183 | |
184 | |
185 | |
186 | |
187 end | 153 end |