Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_compatible_2.m @ 267:f7ac3cd6eeaa operator_remake

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Sparsified all implementation files, removed all matlab warnings, fixed small bugs on minimum grid points.
author Martin Almquist <martin.almquist@it.uu.se>
date Fri, 09 Sep 2016 14:53:41 +0200
parents bfa130b7abf6
children
comparison
equal deleted inserted replaced
266:bfa130b7abf6 267:f7ac3cd6eeaa
24 % randsplutningar f?r D1 och D2 24 % randsplutningar f?r D1 och D2
25 25
26 % Vi b?rjar med normen. Notera att alla SBP operatorer delar samma norm, 26 % Vi b?rjar med normen. Notera att alla SBP operatorer delar samma norm,
27 % vilket ?r n?dv?ndigt f?r stabilitet 27 % vilket ?r n?dv?ndigt f?r stabilitet
28 28
29 BP = 1; 29 BP = 4;
30 if(m<2*BP) 30 if(m<2*BP)
31 error(['Operator requires at least ' num2str(2*BP) ' grid points']); 31 error(['Operator requires at least ' num2str(2*BP) ' grid points']);
32 end 32 end
33 33
34 H=diag(ones(m,1),0);H(1,1)=1/2;H(m,m)=1/2; 34 H=speye(m,m);H(1,1)=1/2;H(m,m)=1/2;
35 35
36 36
37 H=H*h; 37 H=H*h;
38 HI=inv(H); 38 HI=inv(H);
39 39
40 40
41 % First derivative SBP operator, 1st order accurate at first 6 boundary points 41 % First derivative SBP operator, 1st order accurate at first 6 boundary points
42 42
43 q1=1/2; 43 q1=1/2;
44 Q=q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); 44 % Q=q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1));
45 stencil = [-q1,0,q1];
46 d = (length(stencil)-1)/2;
47 diags = -d:d;
48 Q = stripeMatrix(stencil, diags, m);
45 49
46 %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)); 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));
47 51
48 52
49 e_1=zeros(m,1);e_1(1)=1; 53 e_1=sparse(m,1);e_1(1)=1;
50 e_m=zeros(m,1);e_m(m)=1; 54 e_m=sparse(m,1);e_m(m)=1;
51 55
52 56
53 D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ; 57 D1=H\(Q-1/2*(e_1*e_1')+1/2*(e_m*e_m')) ;
54 58
55 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 59 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
56 60
57 61
58 62
59 % Second derivative, 1st order accurate at first 6 boundary points 63 % Second derivative, 1st order accurate at first 6 boundary points
60 m1=-1;m0=2; 64 % m1=-1;m0=2;
61 M=m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0);M(1,1)=1;M(m,m)=1; 65 % % M=m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0);M(1,1)=1;M(m,m)=1;
62 M=M/h; 66 % stencil = [m2,m1,m0,m1,m2];
67 % d = (length(stencil)-1)/2;
68 % diags = -d:d;
69 % M = stripeMatrix(stencil, diags, m);
70 % M=M/h;
63 71
64 S_U=[-1 1]/h; 72 S_U=[-1 1]/h;
65 S_1=zeros(1,m); 73 S_1=sparse(1,m);
66 S_1(1:2)=S_U; 74 S_1(1:2)=S_U;
67 S_m=zeros(1,m); 75 S_m=sparse(1,m);
68 76
69 S_m(m-1:m)=fliplr(-S_U); 77 S_m(m-1:m)=fliplr(-S_U);
70 78
71 D2=HI*(-M-e_1*S_1+e_m*S_m); 79 % D2=H\(-M-e_1*S_1+e_m*S_m);
72 80
73 81
74 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 82 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
75 83
76 84
77 85
78 % Third derivative, 1st order accurate at first 6 boundary points 86 % Third derivative, 1st order accurate at first 6 boundary points
79 87
80 q2=1/2;q1=-1; 88 % q2=1/2;q1=-1;
81 Q3=q2*(diag(ones(m-2,1),2)-diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1)); 89 % % Q3=q2*(diag(ones(m-2,1),2)-diag(ones(m-2,1),-2))+q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1));
90 % stencil = [-q2,-q1,0,q1,q2];
91 % d = (length(stencil)-1)/2;
92 % diags = -d:d;
93 % Q3 = stripeMatrix(stencil, diags, m);
82 94
83 %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)); 95 %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));
84 96
85 97
86 Q3_U = [0 -0.2e1 / 0.5e1 0.3e1 / 0.10e2 0.1e1 / 0.10e2; 0.2e1 / 0.5e1 0 -0.7e1 / 0.10e2 0.3e1 / 0.10e2; -0.3e1 / 0.10e2 0.7e1 / 0.10e2 0 -0.9e1 / 0.10e2; -0.1e1 / 0.10e2 -0.3e1 / 0.10e2 0.9e1 / 0.10e2 0;]; 98 % Q3_U = [0 -0.2e1 / 0.5e1 0.3e1 / 0.10e2 0.1e1 / 0.10e2; 0.2e1 / 0.5e1 0 -0.7e1 / 0.10e2 0.3e1 / 0.10e2; -0.3e1 / 0.10e2 0.7e1 / 0.10e2 0 -0.9e1 / 0.10e2; -0.1e1 / 0.10e2 -0.3e1 / 0.10e2 0.9e1 / 0.10e2 0;];
87 Q3(1:4,1:4)=Q3_U; 99 % Q3(1:4,1:4)=Q3_U;
88 Q3(m-3:m,m-3:m)=flipud( fliplr( -Q3_U ) ); 100 % Q3(m-3:m,m-3:m)=rot90( -Q3_U ,2 );
89 Q3=Q3/h^2; 101 % Q3=Q3/h^2;
90 102
91 103
92 104
93 S2_U=[1 -2 1;]/h^2; 105 S2_U=[1 -2 1;]/h^2;
94 S2_1=zeros(1,m); 106 S2_1=sparse(1,m);
95 S2_1(1:3)=S2_U; 107 S2_1(1:3)=S2_U;
96 S2_m=zeros(1,m); 108 S2_m=sparse(1,m);
97 S2_m(m-2:m)=fliplr(S2_U); 109 S2_m(m-2:m)=fliplr(S2_U);
98 110
99 111
100 112
101 D3=HI*(Q3 - e_1*S2_1 + e_m*S2_m +1/2*S_1'*S_1 -1/2*S_m'*S_m ) ; 113 % D3=H\(Q3 - e_1*S2_1 + e_m*S2_m +1/2*(S_1'*S_1) -1/2*(S_m'*S_m) ) ;
102 114
103 % Fourth derivative, 0th order accurate at first 6 boundary points (still 115 % Fourth derivative, 0th order accurate at first 6 boundary points (still
104 % yield 4th order convergence if stable: for example u_tt=-u_xxxx 116 % yield 4th order convergence if stable: for example u_tt=-u_xxxx
105 117
106 m2=1;m1=-4;m0=6; 118 m2=1;m1=-4;m0=6;
110 122
111 M4_U=[0.4e1 / 0.5e1 -0.7e1 / 0.5e1 0.2e1 / 0.5e1 0.1e1 / 0.5e1; -0.7e1 / 0.5e1 0.16e2 / 0.5e1 -0.11e2 / 0.5e1 0.2e1 / 0.5e1; 0.2e1 / 0.5e1 -0.11e2 / 0.5e1 0.21e2 / 0.5e1 -0.17e2 / 0.5e1; 0.1e1 / 0.5e1 0.2e1 / 0.5e1 -0.17e2 / 0.5e1 0.29e2 / 0.5e1;]; 123 M4_U=[0.4e1 / 0.5e1 -0.7e1 / 0.5e1 0.2e1 / 0.5e1 0.1e1 / 0.5e1; -0.7e1 / 0.5e1 0.16e2 / 0.5e1 -0.11e2 / 0.5e1 0.2e1 / 0.5e1; 0.2e1 / 0.5e1 -0.11e2 / 0.5e1 0.21e2 / 0.5e1 -0.17e2 / 0.5e1; 0.1e1 / 0.5e1 0.2e1 / 0.5e1 -0.17e2 / 0.5e1 0.29e2 / 0.5e1;];
112 124
113 M4(1:4,1:4)=M4_U; 125 M4(1:4,1:4)=M4_U;
114 126
115 M4(m-3:m,m-3:m)=flipud( fliplr( M4_U ) ); 127 M4(m-3:m,m-3:m)=rot90( M4_U ,2 );
116 M4=M4/h^3; 128 M4=M4/h^3;
117 129
118 S3_U=[-1 3 -3 1;]/h^3; 130 S3_U=[-1 3 -3 1;]/h^3;
119 S3_1=zeros(1,m); 131 S3_1=sparse(1,m);
120 S3_1(1:4)=S3_U; 132 S3_1(1:4)=S3_U;
121 S3_m=zeros(1,m); 133 S3_m=sparse(1,m);
122 S3_m(m-3:m)=fliplr(-S3_U); 134 S3_m(m-3:m)=fliplr(-S3_U);
123 135
124 D4=HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m); 136 D4=H\(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m);
125 137
126 138
127 139
128 S_1 = S_1'; 140 S_1 = S_1';
129 S_m = S_m'; 141 S_m = S_m';