-
Notifications
You must be signed in to change notification settings - Fork 0
/
SBP_I8.m
94 lines (67 loc) · 7.37 KB
/
SBP_I8.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
function [S,MM,BD,QQ,H,x,h] = SBP_I8(m,L)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 8:te Implicit SBP Finita differens %%%
%%% operatorer framtagna av Ken Mattsson %%%
%%% %%%
%%% Datum: 2018 06 18 %%%
%%% %%%
%%% 6 randpunkter, bandad norm %%
%%% %%%
%%% H (Normen) %%%
%%% D1 (approx färsta derivatan) %%%
%%% D2 (approx andra derivatan) %%%
%%% S Artificiell Dissipation %%%
%%% %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% With 6 boundary points it was possible to have
% D1 with 5-8-5 accuracy, that had 3 free parameters
% optmised for two leading errors
%
% D2 with 7 boundary points we optmised for 2 leading errors
%L=1;m=21;h=L/(m-1);
% D1 har noggrannhet 5-8-5
% D2 har noggrannhet 3-8-3
%x=linspace(0,L,m);
h=L/(m-1);
x=linspace(0,L,m);
H_U=[0.586619755492542413e18 / 0.3171452440429324800e19 0.441108576863311949e18 / 0.1761918022460736000e19 -0.152467483094515753e18 / 0.660719258422776000e18 0.9046642093687183e16 / 0.49553944381708200e17 -0.165991071774429997e18 / 0.1761918022460736000e19 0.138731627502637093e18 / 0.5285754067382208000e19; 0.441108576863311949e18 / 0.1761918022460736000e19 0.7410635203907926091e19 / 0.5285754067382208000e19 -0.155438554083168289e18 / 0.330359629211388000e18 0.169132174262748227e18 / 0.660719258422776000e18 -0.157501775667810059e18 / 0.5285754067382208000e19 -0.178171110770297671e18 / 0.5285754067382208000e19; -0.152467483094515753e18 / 0.660719258422776000e18 -0.155438554083168289e18 / 0.330359629211388000e18 0.3251789201036087461e19 / 0.1982157775268328000e19 -0.267677573313563593e18 / 0.660719258422776000e18 0.52981923674934827e17 / 0.660719258422776000e18 0.281060957537359e15 / 0.6194243047713525e16; 0.9046642093687183e16 / 0.49553944381708200e17 0.169132174262748227e18 / 0.660719258422776000e18 -0.267677573313563593e18 / 0.660719258422776000e18 0.2469064111296089861e19 / 0.1982157775268328000e19 -0.2712773118379789e16 / 0.330359629211388000e18 -0.15333596630684011e17 / 0.220239752807592000e18; -0.165991071774429997e18 / 0.1761918022460736000e19 -0.157501775667810059e18 / 0.5285754067382208000e19 0.52981923674934827e17 / 0.660719258422776000e18 -0.2712773118379789e16 / 0.330359629211388000e18 0.4738392869842476491e19 / 0.5285754067382208000e19 0.629325798409568167e18 / 0.5285754067382208000e19; 0.138731627502637093e18 / 0.5285754067382208000e19 -0.178171110770297671e18 / 0.5285754067382208000e19 0.281060957537359e15 / 0.6194243047713525e16 -0.15333596630684011e17 / 0.220239752807592000e18 0.629325798409568167e18 / 0.5285754067382208000e19 0.2701816324933024717e19 / 0.3171452440429324800e19;];
h3=1/140;h2=-3/70;h1=3/28;h0=6/7;
H=h3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+h2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+h1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+h0*diag(ones(m,1),0);
H(1:6,1:6)=H_U;
H(m-5:m,m-5:m)=rot90(H_U,2);
H=H*h;
% First derivative SBP operator, 4th order accurate at first 5 boundary points
q3=1/60;q2=-3/20;q1=3/4;
Q=q3*(diag(ones(m-3,1),3) - diag(ones(m-3,1),-3))+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));
Q_U = [0 0.2173462764501499e16 / 0.2473060231152000e16 -0.130333876643261e15 / 0.183189646752000e15 0.850653774171011e15 / 0.1648706820768000e16 -0.163985717405831e15 / 0.706588637472000e15 0.80362690281611e14 / 0.1648706820768000e16; -0.2173462764501499e16 / 0.2473060231152000e16 0 0.303817538718937e15 / 0.197844818492160e15 -0.150024039882917e15 / 0.141317727494400e15 0.19647647570483e14 / 0.36637929350400e14 -0.650103965083537e15 / 0.4946120462304000e16; 0.130333876643261e15 / 0.183189646752000e15 -0.303817538718937e15 / 0.197844818492160e15 0 0.52482306327439e14 / 0.41217670519200e14 -0.614395002823459e15 / 0.989224092460800e15 0.40501050163333e14 / 0.235529545824000e15; -0.850653774171011e15 / 0.1648706820768000e16 0.150024039882917e15 / 0.141317727494400e15 -0.52482306327439e14 / 0.41217670519200e14 0 0.38046246883237e14 / 0.39568963698432e14 -0.413073171846469e15 / 0.1648706820768000e16; 0.163985717405831e15 / 0.706588637472000e15 -0.19647647570483e14 / 0.36637929350400e14 0.614395002823459e15 / 0.989224092460800e15 -0.38046246883237e14 / 0.39568963698432e14 0 0.274844355340637e15 / 0.353294318736000e15; -0.80362690281611e14 / 0.1648706820768000e16 0.650103965083537e15 / 0.4946120462304000e16 -0.40501050163333e14 / 0.235529545824000e15 0.413073171846469e15 / 0.1648706820768000e16 -0.274844355340637e15 / 0.353294318736000e15 0;];
Q(1:6,1:6)=Q_U;
Q(m-5:m,m-5:m)=rot90(-Q_U,2);
e_l=zeros(m,1);e_l(1)=1;
e_r=zeros(m,1);e_r(m)=1;
QQ=Q-1/2*e_l*e_l'+1/2*e_r*e_r';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Second derivative, 2nd order accurate at first 6 boundary points
m4=-3/560;m3=2/63;m2=0;m1=-6/5;m0=169/72;
M=m4*(diag(ones(m-4,1),4)+diag(ones(m-4,1),-4))+m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0);
% Below with 6 boundary points There was one free parameter left,
% optimised for leading order arror
M_U=[0.13770968800936061e17 / 0.9892240924608000e16 -0.4300348585643179e16 / 0.1978448184921600e16 0.66937433609903e14 / 0.47105909164800e14 -0.951197030290193e15 / 0.989224092460800e15 0.779545212871249e15 / 0.1978448184921600e16 -0.79093632472679e14 / 0.1099137880512000e16; -0.4300348585643179e16 / 0.1978448184921600e16 0.1270292308260497e16 / 0.219827576102400e15 -0.5870518615787041e16 / 0.989224092460800e15 0.3455893930366859e16 / 0.989224092460800e15 -0.943081398635729e15 / 0.659482728307200e15 0.526211378046257e15 / 0.1978448184921600e16; 0.66937433609903e14 / 0.47105909164800e14 -0.5870518615787041e16 / 0.989224092460800e15 0.4464641360686277e16 / 0.494612046230400e15 -0.116090342903909e15 / 0.18318964675200e14 0.2198253388831819e16 / 0.989224092460800e15 -0.388525668633169e15 / 0.989224092460800e15; -0.951197030290193e15 / 0.989224092460800e15 0.3455893930366859e16 / 0.989224092460800e15 -0.116090342903909e15 / 0.18318964675200e14 0.3186195709270117e16 / 0.494612046230400e15 -0.2984932746782561e16 / 0.989224092460800e15 0.116872806771529e15 / 0.329741364153600e15; 0.779545212871249e15 / 0.1978448184921600e16 -0.943081398635729e15 / 0.659482728307200e15 0.2198253388831819e16 / 0.989224092460800e15 -0.2984932746782561e16 / 0.989224092460800e15 0.2087760214166131e16 / 0.659482728307200e15 -0.2692431992885291e16 / 0.1978448184921600e16; -0.79093632472679e14 / 0.1099137880512000e16 0.526211378046257e15 / 0.1978448184921600e16 -0.388525668633169e15 / 0.989224092460800e15 0.116872806771529e15 / 0.329741364153600e15 -0.2692431992885291e16 / 0.1978448184921600e16 0.23531662112543101e17 / 0.9892240924608000e16;];
M(1:6,1:6)=M_U;
M(m-5:m,m-5:m)=rot90(M_U,2);
M=M/h;
d1_U=[-137/60 5 -5 10/3 -5/4 1/5]/h;
d1_l=zeros(1,m);
d1_l(1:6)=d1_U;
d1_r=zeros(1,m);
d1_r(m-5:m)=fliplr(-d1_U);
BD=-e_l*d1_l+e_r*d1_r;
MM=-M+BD;
% Artificial dissipation = -HI*S (in RHS)
% and will lead to order 4-9-4
d5=0*[-1 5 -10 10 -5 1];
DD_5=(-diag(ones(m-3,1),-3)+5*diag(ones(m-2,1),-2)-10*diag(ones(m-1,1),-1)+10*diag(ones(m,1),0)-5*diag(ones(m-1,1),1)+diag(ones(m-2,1),2));
DD_5(1:3,1:6)=[d5;d5;d5];
DD_5(m-1:m,m-5:m)=[d5;d5];
DD=DD_5'*DD_5;
aa=max(max(DD));
S=-1/aa*DD;