Mercurial > repos > public > sbplib

comparison +sbp/higher_variable4.m @ 249:02423f9323c6 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Started modifying operator files.
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 05 Sep 2016 16:38:41 +0200
parents c2ca9717db4d
children
comparison
equal deleted inserted replaced
248:81e0ead29431 249:02423f9323c6
5 %%% H (Normen) %%% 5 %%% H (Normen) %%%
6 %%% D1=H^(-1)Q (approx f?rsta derivatan) %%% 6 %%% D1=H^(-1)Q (approx f?rsta derivatan) %%%
7 %%% D2 (approx andra derivatan) %%% 7 %%% D2 (approx andra derivatan) %%%
8 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 8 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
9 9
10 Hv = ones(m,1);
11 Hv(1:4) = [17/48 59/48 43/48 49/48];
12 Hv(m-3:m) = rot90(Hv(1:4),2);
13 Hv = h * Hv;
14 H = spdiag(Hv,0);
15 HI = spdiag(1./Hv,0);
10 16
11 %m=20; %problemstorlek 17 e_1 = sparse(m, 1);
12 %h=1/(m-1); 18 e_1(1) = 1;
13 %h=1; 19 e_m = rot90(e_1, 2);
14 20
15 c = ones(m,1); 21 S_1 = sparse(m, 1);
22 S_1(1:4) = [-0.11e2/0.6e1 3 -0.3e1/0.2e1 0.1e1/0.3e1]/h;
23 S_m = -rot90(S_1, 2);
24
25 S2_1 = sparse(m, 1);
26 S2_1(1:4) = [2 -5 4 -1;]/h^2;
27 S2_m = rot90(S2_1, 2);
28
29 S3_1 = sparse(m, 1);
30 S3_1(1:4) = [-1 3 -3 1;]/h^3;
31 S3_m = rot90(S3_1, 2);
16 32
17 33
18 H=diag(ones(m,1),0); 34 q_diags = [-2, -1, 0, 1, 2];
19 H(1:4,1:4)=diag([17/48 59/48 43/48 49/48]); 35 q_stencil = [-1/12, 8/12, 0, -8/12, 1/12];
20 H(m-3:m,m-3:m)=fliplr(flipud(diag([17/48 59/48 43/48 49/48]))); 36 Q = stripeMatrix(q_stencil, q_diags,m);
21 H=H*h;
22 HI=inv(H);
23 HI = sparse(HI);
24 37
38 Q_U = [
39 0 0.59e2/0.96e2 -0.1e1/0.12e2 -0.1e1/0.32e2;
40 -0.59e2/0.96e2 0 0.59e2/0.96e2 0;
41 0.1e1/0.12e2 -0.59e2/0.96e2 0 0.59e2/0.96e2;
42 0.1e1/0.32e2 0 -0.59e2/0.96e2 0;
43 ];
25 44
45 Q(1:4,1:4)=Q_U;
46 Q(m-3:m,m-3:m)=-rot90(Q_U, 2);
26 47
27 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));
28 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;];
29 Q(1:4,1:4)=Q_U;
30 Q(m-3:m,m-3:m)=flipud( fliplr(-Q_U(1:4,1:4) ) );
31
32 e_1=zeros(m,1);e_1(1)=1;
33 e_m=zeros(m,1);e_m(m)=1;
34 48
35 D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ; 49 D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ;
36 50
37 M_U=[0.9e1 / 0.8e1 -0.59e2 / 0.48e2 0.1e1 / 0.12e2 0.1e1 / 0.48e2; -0.59e2 / 0.48e2 0.59e2 / 0.24e2 -0.59e2 / 0.48e2 0; 0.1e1 / 0.12e2 -0.59e2 / 0.48e2 0.55e2 / 0.24e2 -0.59e2 / 0.48e2; 0.1e1 / 0.48e2 0 -0.59e2 / 0.48e2 0.59e2 / 0.24e2;]; 51 m_diags = [-2, -1, 0, 1, 2];
38 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)); 52 m_stencil = [-1/12, 16/12, -30/12, -16/12, -1/12];
53
54 M_U=[
55 0.9e1 / 0.8e1 -0.59e2 / 0.48e2 0.1e1 / 0.12e2 0.1e1 / 0.48e2;
56 -0.59e2 / 0.48e2 0.59e2 / 0.24e2 -0.59e2 / 0.48e2 0;
57 0.1e1 / 0.12e2 -0.59e2 / 0.48e2 0.55e2 / 0.24e2 -0.59e2 / 0.48e2;
58 0.1e1 / 0.48e2 0 -0.59e2 / 0.48e2 0.59e2 / 0.24e2;
59 ];
39 60
40 M(1:4,1:4)=M_U; 61 M(1:4,1:4)=M_U;
41 62 M(m-3:m,m-3:m)=rot90(M_U, 2);
42 M(m-3:m,m-3:m)=flipud( fliplr( M_U ) );
43 M=M/h; 63 M=M/h;
44
45 S_U=[-0.11e2 / 0.6e1 3 -0.3e1 / 0.2e1 0.1e1 / 0.3e1;]/h;
46 S_1=zeros(1,m);
47 S_1(1:4)=S_U;
48 S_m=zeros(1,m);
49 S_m(m-3:m)=fliplr(-S_U);
50 S_1 = S_1';
51 S_m = S_m';
52
53
54 M=sparse(m,m);
55 e_1 = sparse(e_1);
56 e_m = sparse(e_m);
57 S_1 = sparse(S_1);
58 S_m = sparse(S_m);
59 64
60 65
61 scheme_width = 5; 66 scheme_width = 5;
62 scheme_radius = (scheme_width-1)/2; 67 scheme_radius = (scheme_width-1)/2;
63 r = (1+scheme_radius):(m-scheme_radius); 68 r = (1+scheme_radius):(m-scheme_radius);
114 M=M/h; 119 M=M/h;
115 D2=HI*(-M-c(1)*e_1*S_1'+c(m)*e_m*S_m'); 120 D2=HI*(-M-c(1)*e_1*S_1'+c(m)*e_m*S_m');
116 end 121 end
117 D2 = @D2_fun; 122 D2 = @D2_fun;
118 123
119
120 S2_U=[2 -5 4 -1;]/h^2;
121 S2_1=zeros(1,m);
122 S2_1(1:4)=S2_U;
123 S2_m=zeros(1,m);
124 S2_m(m-3:m)=fliplr(S2_U);
125 S2_1 = S2_1';
126 S2_m = S2_m';
127
128 m3=-1/6;m2=2;m1=-13/2;m0=28/3; 124 m3=-1/6;m2=2;m1=-13/2;m0=28/3;
129 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); 125 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);
130 126
131 %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)); 127 %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));
132 128
135 M4(1:6,1:6)=M4_U; 131 M4(1:6,1:6)=M4_U;
136 132
137 M4(m-5:m,m-5:m)=flipud( fliplr( M4_U ) ); 133 M4(m-5:m,m-5:m)=flipud( fliplr( M4_U ) );
138 M4=M4/h^3; 134 M4=M4/h^3;
139 135
140 S3_U=[-1 3 -3 1;]/h^3;
141 S3_1=zeros(1,m);
142 S3_1(1:4)=S3_U;
143 S3_m=zeros(1,m);
144 S3_m(m-3:m)=fliplr(-S3_U);
145 S3_1 = S3_1';
146 S3_m = S3_m';
147 136
148 D4=HI*(M4-e_1*S3_1'+e_m*S3_m' + S_1*S2_1'-S_m*S2_m'); 137 D4=HI*(M4-e_1*S3_1'+e_m*S3_m' + S_1*S2_1'-S_m*S2_m');
149 138
150 139
151 140