Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_variable_2.m @ 312:9230c056a574 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fixed formatting.
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 23 Sep 2016 19:14:04 +0200
parents 713b125038a3
children 52b4cdf27633
comparison
equal deleted inserted replaced
311:713b125038a3 312:9230c056a574
1 % Returns D2 as a function handle 1 % Returns D2 as a function handle
2 function [H, HI, D1, D2, D3, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = d4_variable_2(m,h) 2 function [H, HI, D1, D2, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = d4_variable_2(m,h)
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 %%% 4:de ordn. SBP Finita differens %%% 4 %%% 4:de ordn. SBP Finita differens %%%
5 %%% operatorer framtagna av Ken Mattsson %%% 5 %%% operatorer framtagna av Ken Mattsson %%%
6 %%% %%% 6 %%% %%%
7 %%% 6 randpunkter, diagonal norm %%% 7 %%% 6 randpunkter, diagonal norm %%%
8 %%% %%% 8 %%% %%%
9 %%% Datum: 2013-11-11 %%% 9 %%% Datum: 2013-11-11 %%%
10 %%% %%%
11 %%% %%%
12 %%% H (Normen) %%%
13 %%% D1 (approx f?rsta derivatan) %%%
14 %%% D2 (approx andra derivatan) %%%
15 %%% D3 (approx tredje derivatan) %%%
16 %%% D2 (approx fj?rde derivatan) %%%
17 %%% %%%
18 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 10 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
19
20 % M?ste ange antal punkter (m) och stegl?ngd (h)
21 % Notera att dessa opetratorer ?r framtagna f?r att anv?ndas n?r
22 % vi har 3de och 4de derivator i v?r PDE
23 % I annat fall anv?nd de "traditionella" som har noggrannare
24 % randsplutningar f?r D1 och D2
25
26 % Vi b?rjar med normen. Notera att alla SBP operatorer delar samma norm,
27 % vilket ?r n?dv?ndigt f?r stabilitet
28 11
29 BP = 4; 12 BP = 4;
30 if(m<2*BP) 13 if(m<2*BP)
31 error(['Operator requires at least ' num2str(2*BP) ' grid points']); 14 error(['Operator requires at least ' num2str(2*BP) ' grid points']);
32 end 15 end
33 16
34 H=speye(m,m);H(1,1)=1/2;H(m,m)=1/2;
35 17
18 H=speye(m,m);
19 H(1,1)=1/2;
20 H(m,m)=1/2;
36 21
37 H=H*h; 22 H=H*h;
38 HI=inv(H); 23 HI=inv(H);
39 24
40 25
47 diags = -d:d; 32 diags = -d:d;
48 Q = stripeMatrix(stencil, diags, m); 33 Q = stripeMatrix(stencil, diags, m);
49 34
50 %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)); 35 %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));
51 36
52 37 e_1=sparse(m,1);
53 e_1=sparse(m,1);e_1(1)=1; 38 e_1(1)=1;
54 e_m=sparse(m,1);e_m(m)=1; 39 e_m=sparse(m,1);
55 40 e_m(m)=1;
56 41
57 D1=HI*(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ; 42 D1=HI*(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ;
58
59 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
60
61 43
62 44
63 % Second derivative, 1st order accurate at first boundary points 45 % Second derivative, 1st order accurate at first boundary points
64 46
65 % below for constant coefficients 47 % below for constant coefficients
125 Q3 = stripeMatrix(stencil, diags, m); 107 Q3 = stripeMatrix(stencil, diags, m);
126 108
127 %QQ3=(-1/8*diag(ones(m-3,1),3) + 1*diag(ones(m-2,1),2) - 13/8*diag(ones(m-1,1),1) +13/8*diag(ones(m-1,1),-1) -1*diag(ones(m-2,1),-2) + 1/8*diag(ones(m-3,1),-3)); 109 %QQ3=(-1/8*diag(ones(m-3,1),3) + 1*diag(ones(m-2,1),2) - 13/8*diag(ones(m-1,1),1) +13/8*diag(ones(m-1,1),-1) -1*diag(ones(m-2,1),-2) + 1/8*diag(ones(m-3,1),-3));
128 110
129 111
130 Q3_U = [0 -0.13e2 / 0.16e2 0.7e1 / 0.8e1 -0.1e1 / 0.16e2; 0.13e2 / 0.16e2 0 -0.23e2 / 0.16e2 0.5e1 / 0.8e1; -0.7e1 / 0.8e1 0.23e2 / 0.16e2 0 -0.17e2 / 0.16e2; 0.1e1 / 0.16e2 -0.5e1 / 0.8e1 0.17e2 / 0.16e2 0;]; 112 Q3_U = [
113 0 -0.13e2/0.16e2 0.7e1/0.8e1 -0.1e1/0.16e2;
114 0.13e2/0.16e2 0 -0.23e2/0.16e2 0.5e1/0.8e1;
115 -0.7e1/0.8e1 0.23e2/0.16e2 0 -0.17e2/0.16e2;
116 0.1e1/0.16e2 -0.5e1/0.8e1 0.17e2/0.16e2 0;
117 ];
131 Q3(1:4,1:4)=Q3_U; 118 Q3(1:4,1:4)=Q3_U;
132 Q3(m-3:m,m-3:m)=rot90( -Q3_U ,2 ); 119 Q3(m-3:m,m-3:m)=rot90( -Q3_U ,2 );
133 Q3=Q3/h^2; 120 Q3=Q3/h^2;
134 121
135 122
156 diags = -d:d; 143 diags = -d:d;
157 M4 = stripeMatrix(stencil, diags, m); 144 M4 = stripeMatrix(stencil, diags, m);
158 145
159 %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)); 146 %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));
160 147
161 M4_U=[0.13e2 / 0.10e2 -0.12e2 / 0.5e1 0.9e1 / 0.10e2 0.1e1 / 0.5e1; -0.12e2 / 0.5e1 0.26e2 / 0.5e1 -0.16e2 / 0.5e1 0.2e1 / 0.5e1; 0.9e1 / 0.10e2 -0.16e2 / 0.5e1 0.47e2 / 0.10e2 -0.17e2 / 0.5e1; 0.1e1 / 0.5e1 0.2e1 / 0.5e1 -0.17e2 / 0.5e1 0.29e2 / 0.5e1;]; 148 M4_U=[
162 149 0.13e2/0.10e2 -0.12e2/0.5e1 0.9e1/0.10e2 0.1e1/0.5e1;
150 -0.12e2/0.5e1 0.26e2/0.5e1 -0.16e2/0.5e1 0.2e1/0.5e1;
151 0.9e1/0.10e2 -0.16e2/0.5e1 0.47e2/0.10e2 -0.17e2/0.5e1;
152 0.1e1/0.5e1 0.2e1/0.5e1 -0.17e2/0.5e1 0.29e2/0.5e1;
153 ];
163 154
164 M4(1:4,1:4)=M4_U; 155 M4(1:4,1:4)=M4_U;
165 156
166 M4(m-3:m,m-3:m)=rot90( M4_U ,2 ); 157 M4(m-3:m,m-3:m)=rot90( M4_U ,2 );
167 M4=M4/h^3; 158 M4=M4/h^3;