Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d1_gauss_4.m @ 405:4d9d8064e58b feature/SBPInTimeGauss

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implementation of D1 based on Gauss quadrature formula with 4 nodes.
author Martin Almquist <martin.almquist@it.uu.se>
date Thu, 02 Feb 2017 17:05:43 +0100
parents
children ba73c9c8d1a6
comparison
equal deleted inserted replaced
404:d6d27fdc342a 405:4d9d8064e58b
1 function [D1,H,x,h,e_l,e_r] = d1_gauss_4(N,L)
2
3 % L: Domain length
4 % N: Number of grid points
5 if(nargin < 2)
6 L = 1;
7 end
8
9 if(N~=4)
10 error('This operator requires exactly 4 grid points');
11 end
12
13 % Quadrature nodes on interval [-1, 1]
14 x = [ -0.8611363115940526; -0.3399810435848563; 0.3399810435848563; 0.8611363115940526];
15
16 % Shift nodes to [0,L]
17 x = (x+1)/2*L;
18
19 % Boundary extrapolation operators
20 e_l = [1.5267881254572668; -0.8136324494869273; 0.4007615203116504; -0.1139171962819899];
21 e_r = flipud(e_l);
22 e_l = sparse(e_l);
23 e_r = sparse(e_r);
24
25 %%%% Compute approximate h %%%%%%%%%%
26 h = L/(N-1);
27 %%%%%%%%%%%%%%%%%%%%%%%%%
28
29 %%%% Norm matrix on [-1,1] %%%%%%%%
30 P = sparse(N,N);
31 P(1,1) = 0.3478548451374539;
32 P(2,2) = 0.6521451548625461;
33 P(3,3) = 0.6521451548625461;
34 P(4,4) = 0.3478548451374539;
35 %%%%%%%%%%%%%%%%%%%%%%%%%
36
37 %%%% Norm matrix on [0,L] %%%%%%%%
38 H = P*L/2;
39 %%%%%%%%%%%%%%%%%%%%%%%%%
40
41 %%%% D1 on [-1,1] %%%%%%%%
42 D1 = sparse(N,N);
43 D1(1,1) = -3.3320002363522817;
44 D1(1,2) = 4.8601544156851962;
45 D1(1,3) = -2.1087823484951789;
46 D1(1,4) = 0.5806281691622644;
47
48 D1(2,1) = -0.7575576147992339;
49 D1(2,2) = -0.3844143922232086;
50 D1(2,3) = 1.4706702312807167;
51 D1(2,4) = -0.3286982242582743;
52
53 D1(3,1) = 0.3286982242582743;
54 D1(3,2) = -1.4706702312807167;
55 D1(3,3) = 0.3844143922232086;
56 D1(3,4) = 0.7575576147992339;
57
58 D1(4,1) = -0.5806281691622644;
59 D1(4,2) = 2.1087823484951789;
60 D1(4,3) = -4.8601544156851962;
61 D1(4,4) = 3.3320002363522817;
62 %%%%%%%%%%%%%%%%%%%%%%%%%
63
64 % D1 on [0,L]
65 D1 = D1*2/L;