Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d2_variable_2.m @ 362:ded4156e53e2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added 2nd order accurate 2nd derivative with variable coefficents in a separate implementation file, used by the class D2Variable.
author Martin Almquist <martin.almquist@it.uu.se>
date Tue, 20 Dec 2016 15:04:02 +0100
parents
children b758d1cf4c8e
comparison
equal deleted inserted replaced
303:f18142c1530b 362:ded4156e53e2
1 function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_2(m,h)
2
3 BP = 1;
4 if(m<2*BP)
5 error(['Operator requires at least ' num2str(2*BP) ' grid points']);
6 end
7
8 % Norm
9 Hv = ones(m,1);
10 Hv(1) = 1/2;
11 Hv(m:m) = 1/2;
12 Hv = h*Hv;
13 H = spdiag(Hv, 0);
14 HI = spdiag(1./Hv, 0);
15
16
17 % Boundary operators
18 e_l = sparse(m,1);
19 e_l(1) = 1;
20 e_r = rot90(e_l, 2);
21
22 d1_l = sparse(m,1);
23 d1_l(1:3) = 1/h*[-3/2 2 -1/2];
24 d1_r = -rot90(d1_l, 2);
25
26 % D1 operator
27 diags = -1:1;
28 stencil = [-1/2 0 1/2];
29 D1 = stripeMatrix(stencil, diags, m);
30
31 D1(1,1)=-1;D1(1,2)=1;D1(m,m-1)=-1;D1(m,m)=1;
32 D1(m,m-1)=-1;D1(m,m)=1;
33 D1=D1/h;
34 %Q=H*D1 + 1/2*(e_1*e_1') - 1/2*(e_m*e_m');
35
36
37 M=sparse(m,m);
38
39 scheme_width = 3;
40 scheme_radius = (scheme_width-1)/2;
41 r = (1+scheme_radius):(m-scheme_radius);
42
43 function D2 = D2_fun(c)
44
45 Mm1 = -c(r-1)/2 - c(r)/2;
46 M0 = c(r-1)/2 + c(r) + c(r+1)/2;
47 Mp1 = -c(r)/2 - c(r+1)/2;
48
49 M(r,:) = spdiags([Mm1 M0 Mp1],0:2*scheme_radius,length(r),m);
50
51
52 M(1:2,1:2)=[c(1)/2 + c(2)/2 -c(1)/2 - c(2)/2; -c(1)/2 - c(2)/2 c(1)/2 + c(2) + c(3)/2;];
53 M(m-1:m,m-1:m)=[c(m-2)/2 + c(m-1) + c(m)/2 -c(m-1)/2 - c(m)/2; -c(m-1)/2 - c(m)/2 c(m-1)/2 + c(m)/2;];
54 M=M/h;
55
56 D2=HI*(-M-c(1)*e_l*d1_l'+c(m)*e_r*d1_r');
57 end
58 D2 = @D2_fun;
59 end