Mercurial > repos > public > sbplib

changeset 382:8b8672134be8 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Edit for clarity.
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 02 Jan 2017 09:12:05 +0100
parents d32c749408cb
children 151b08366d63
files +sbp/+implementations/d2_2.m +sbp/+implementations/d2_4.m
diffstat 2 files changed, 38 insertions(+), 21 deletions(-) [+]
line wrap: on
line diff
diff -r d32c749408cb -r 8b8672134be8 +sbp/+implementations/d2_2.m
--- a/+sbp/+implementations/d2_2.m	Mon Jan 02 09:07:49 2017 +0100
+++ b/+sbp/+implementations/d2_2.m	Mon Jan 02 09:12:05 2017 +0100
@@ -1,5 +1,5 @@
 function [H, HI, D1, D2, e_1, e_m, M, Q, S_1, S_m] = d2_2(m,h)
-    
+
     BP = 1;
     if(m<2*BP)
         error(['Operator requires at least ' num2str(2*BP) ' grid points']);
@@ -11,33 +11,40 @@
     H=(speye(m,m));H(1,1)=0.5;H(m,m)=0.5;
     H=h*H;
     HI=inv(H);
-    
+
     diags   = -1:1;
     stencil = [-1/2 0 1/2];
     D1 = stripeMatrix(stencil, diags, m);
-    
-    D1(1,1)=-1;D1(1,2)=1;D1(m,m-1)=-1;D1(m,m)=1;
-    D1(m,m-1)=-1;D1(m,m)=1;
-    D1=D1/h;
+
+    D1(1,1) = -1;
+    D1(1,2) = 1;
 
-    Q=H*D1 + 1/2*(e_1*e_1') - 1/2*(e_m*e_m');
+    D1(m,m-1) = -1;
+    D1(m,m)   = 1;
+
+    D1 = D1/h;
+
+    Q = H*D1 + 1/2*(e_1*e_1') - 1/2*(e_m*e_m');
 
     diags   = -1:1;
     stencil = [1 -2 1];
     D2 = stripeMatrix(stencil, diags, m);
-    
-    D2(1,1)=1;D2(1,2)=-2;D2(1,3)=1;
-    D2(m,m-2)=1;D2(m,m-1)=-2;D2(m,m)=1;
-    D2=D2/h^2;
 
-    S_U=[-3/2, 2, -1/2]/h;
-    S_1=sparse(1,m);
-    S_1(1:3)=S_U;
-    S_m=sparse(1,m);
-    S_m(m-2:m)=fliplr(-S_U);
+    D2(1,1) = 1;
+    D2(1,2) = -2;
+    D2(1,3) = 1;
+    D2(m,m-2) = 1;
+    D2(m,m-1) = -2;
+    D2(m,m)   = 1;
+    D2 = D2/h^2;
 
+    S_U = [-3/2, 2, -1/2]/h;
+    S_1 = sparse(1,m);
+    S_1(1:3) = S_U;
+    S_m = sparse(1,m);
+    S_m(m-2:m) = fliplr(-S_U);
 
-    M=-H*D2-e_1*S_1+e_m*S_m;
+    M = -H*D2-e_1*S_1+e_m*S_m;
     S_1 = S_1';
     S_m = S_m';
 end
\ No newline at end of file
diff -r d32c749408cb -r 8b8672134be8 +sbp/+implementations/d2_4.m
--- a/+sbp/+implementations/d2_4.m	Mon Jan 02 09:07:49 2017 +0100
+++ b/+sbp/+implementations/d2_4.m	Mon Jan 02 09:12:05 2017 +0100
@@ -1,5 +1,5 @@
 function [H, HI, D1, D2, e_1, e_m, M, Q, S_1, S_m] = d2_4(m,h)
-    
+
     BP = 4;
     if(m<2*BP)
         error(['Operator requires at least ' num2str(2*BP) ' grid points']);
@@ -20,7 +20,12 @@
 
     Q=spdiags([e -8*e 0*e 8*e -e], -2:2, m, m)/12;
     %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));
-    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;];
+    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;
+    ];
     Q(1:4,1:4)=Q_U;
     Q(m-3:m,m-3:m)=rot90( -Q_U(1:4,1:4) ,2 );
 
@@ -30,13 +35,18 @@
 
     %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));
 
-    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;];
+    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;
+    ];
     M(1:4,1:4)=M_U;
 
     M(m-3:m,m-3:m)=rot90(  M_U ,2 );
     M=M/h;
 
-    S_U=[-0.11e2 / 0.6e1 3 -0.3e1 / 0.2e1 0.1e1 / 0.3e1;]/h;
+    S_U=[-0.11e2/0.6e1 3 -0.3e1/0.2e1 0.1e1/0.3e1;]/h;
     S_1=sparse(sparse(1,m));
     S_1(1:4)=S_U;
     S_m=sparse(sparse(1,m));