Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d2_2.m @ 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 f7ac3cd6eeaa
children
comparison
equal deleted inserted replaced
381:d32c749408cb 382:8b8672134be8
1 function [H, HI, D1, D2, e_1, e_m, M, Q, S_1, S_m] = d2_2(m,h) 1 function [H, HI, D1, D2, e_1, e_m, M, Q, S_1, S_m] = d2_2(m,h)
2 2
3 BP = 1; 3 BP = 1;
4 if(m<2*BP) 4 if(m<2*BP)
5 error(['Operator requires at least ' num2str(2*BP) ' grid points']); 5 error(['Operator requires at least ' num2str(2*BP) ' grid points']);
6 end 6 end
7 7
9 e_m=sparse(m,1);e_m(m)=1; 9 e_m=sparse(m,1);e_m(m)=1;
10 10
11 H=(speye(m,m));H(1,1)=0.5;H(m,m)=0.5; 11 H=(speye(m,m));H(1,1)=0.5;H(m,m)=0.5;
12 H=h*H; 12 H=h*H;
13 HI=inv(H); 13 HI=inv(H);
14 14
15 diags = -1:1; 15 diags = -1:1;
16 stencil = [-1/2 0 1/2]; 16 stencil = [-1/2 0 1/2];
17 D1 = stripeMatrix(stencil, diags, m); 17 D1 = stripeMatrix(stencil, diags, m);
18
19 D1(1,1)=-1;D1(1,2)=1;D1(m,m-1)=-1;D1(m,m)=1;
20 D1(m,m-1)=-1;D1(m,m)=1;
21 D1=D1/h;
22 18
23 Q=H*D1 + 1/2*(e_1*e_1') - 1/2*(e_m*e_m'); 19 D1(1,1) = -1;
20 D1(1,2) = 1;
21
22 D1(m,m-1) = -1;
23 D1(m,m) = 1;
24
25 D1 = D1/h;
26
27 Q = H*D1 + 1/2*(e_1*e_1') - 1/2*(e_m*e_m');
24 28
25 diags = -1:1; 29 diags = -1:1;
26 stencil = [1 -2 1]; 30 stencil = [1 -2 1];
27 D2 = stripeMatrix(stencil, diags, m); 31 D2 = stripeMatrix(stencil, diags, m);
28
29 D2(1,1)=1;D2(1,2)=-2;D2(1,3)=1;
30 D2(m,m-2)=1;D2(m,m-1)=-2;D2(m,m)=1;
31 D2=D2/h^2;
32 32
33 S_U=[-3/2, 2, -1/2]/h; 33 D2(1,1) = 1;
34 S_1=sparse(1,m); 34 D2(1,2) = -2;
35 S_1(1:3)=S_U; 35 D2(1,3) = 1;
36 S_m=sparse(1,m); 36 D2(m,m-2) = 1;
37 S_m(m-2:m)=fliplr(-S_U); 37 D2(m,m-1) = -2;
38 D2(m,m) = 1;
39 D2 = D2/h^2;
38 40
41 S_U = [-3/2, 2, -1/2]/h;
42 S_1 = sparse(1,m);
43 S_1(1:3) = S_U;
44 S_m = sparse(1,m);
45 S_m(m-2:m) = fliplr(-S_U);
39 46
40 M=-H*D2-e_1*S_1+e_m*S_m; 47 M = -H*D2-e_1*S_1+e_m*S_m;
41 S_1 = S_1'; 48 S_1 = S_1';
42 S_m = S_m'; 49 S_m = S_m';
43 end 50 end