-
Notifications
You must be signed in to change notification settings - Fork 0
/
SBP_SL4.m
112 lines (78 loc) · 5.34 KB
/
SBP_SL4.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
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
function [S,MM,BD,QQ,H,x,h] = SBP_SL4(m,L)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 4:de Spectral SBP Finita differens %%%
%%% operatorer framtagna av Ken Mattsson %%%
%%% %%%
%%% Datum: 2019 07 09 %%%
%%% %%%
%%% 4 randpunkter, %%%
%%% %%%
%%% Opt. för w=1.3, 1.7, 2.0, 2.3, 2.5 %%%
%%% %%%
%%% Neptadiagonal norm %%%
%%% Neptadiagonal Q %%%
%%% %%%
%%% Noggrannhet 3-4-3 %%%
%%% %%%
%%% %%%
%%% %%%
%%% H (Normen) %%%
%%% D1 (approx första derivatan) %%%
%%% D2 (approx andra derivatan) %%%
%%% S Artificiell Dissipation %%%
%%% %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% D1 har noggrannhet 3-4-3
%x=linspace(0,L,m);
h=L/(m-1);
x=linspace(0,L,m)';
H_U=[0.272765852144443500436871206195e0 0.110503577576841854695688450492e0 -0.214935939776929923799817276532e-2 -0.123099486212404624101095251785e0; 0.110503577576841854695688450492e0 0.102672568170820815832095910978e1 0.192654317376357480005287217945e0 0.218484672079950056490071930980e0; -0.214935939776929923799817276532e-2 0.192654317376357480005287217945e0 0.234520342145131697556064115881e0 0.210000000000000000000000000000e0; -0.123099486212404624101095251785e0 0.218484672079950056490071930980e0 0.210000000000000000000000000000e0 0.574886914610060495860000000000e0;];
h0 = 0.522989416184906449839043502047;
h1 = 0.271529435853391430058862861221;
h2 = -0.00765473902661188212140935422742;
h3 = -0.0235714284229686377340540977080;
h4 = -0.00179797649626413512292116030866;
H=h4*(diag(ones(m-4,1),4)+diag(ones(m-4,1),-4))+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:4,1:4)=H_U;
H(m-3:m,m-3:m)=rot90(H_U,2);
H=H*h;
q1 = 0.386538021757943806903979515289;
q2 = 0.120984863913399451287779422734;
q3 = -0.0334146219186068540222094702356;
q4 = -0.00706597095723053685322748751251;
Q=q4*(diag(ones(m-4,1),4) - diag(ones(m-4,1),-4))+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.596223946079145138601586516315e0 0.442635292696318892730544080289e-1 -0.133421504391546491021413436832e0; -0.596223946079145138601586516315e0 0 0.218126430582242711421063652503e0 0.418578108372739818055959821561e0; -0.442635292696318892730544080289e-1 -0.218126430582242711421063652503e0 0 0.181885688814312540281775595545e0; 0.133421504391546491021413436832e0 -0.418578108372739818055959821561e0 -0.181885688814312540281775595545e0 0;];
Q(1:4,1:4)=Q_U;
Q(m-3:m,m-3: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 SBP
% Second derivative, 2nd order accurate at first 5 boundary points
% This one optimized dispersion with 2 adional free coefficients , chosen to be exactat w=2 and w=2.6
m0 = 0.636917370013909589036000761283;
m1 = -0.0142217417023750211989443636436;
m2 = -0.319876406618387777350327110712;
m3 = -0.00584636180440348467230133436450;
m4 = 0.0212001115044144138820384317776;
m5 = 0.000285713613797074821533996294469;
M=m5*(diag(ones(m-5,1),5)+diag(ones(m-5,1),-5))+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 5 boundary points There was one free parameter left outher most coeficient in interior stencil,
% optimised for leading order arror
M_U=[0.108429546134659364757697328033e1 -0.104251612712475322873993929885e1 -0.134732707933946570734100853615e0 0.601176654695753122016604787183e-1 0.325499946287337648738723971147e-1; -0.104251612712475322873993929885e1 0.158762367951981935941376440505e1 -0.528303719788691208033292938948e-1 -0.465374247370698358149771007596e0 -0.483887581637101404242972327766e-1; -0.134732707933946570734100853615e0 -0.528303719788691208033292938948e-1 0.305417952970789357566606515851e0 0.123997043776631497769555782530e0 -0.257491380148413167830003244562e0; 0.601176654695753122016604787183e-1 -0.465374247370698358149771007596e0 0.123997043776631497769555782530e0 0.640832051421767624340881043927e0 -0.553355699926963028432702805769e-1; 0.325499946287337648738723971147e-1 -0.483887581637101404242972327766e-1 -0.257491380148413167830003244562e0 -0.553355699926963028432702805769e-1 0.647124398683040640741698741437e0;];
M(1:5,1:5)=M_U;
M(m-4:m,m-4:m)=rot90(M_U,2);
M=M/h;
d1_U=[-11/6,3, -3/2, 1/3]/h;
d1_l=zeros(1,m);
d1_l(1:4)=d1_U;
d1_r=zeros(1,m);
d1_r(m-3: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 3-7-3
d4=0*[1 -4 6 -4 1];
DD_4=(diag(ones(m-2,1),2)-4*diag(ones(m-1,1),1)+6*diag(ones(m,1),0)-4*diag(ones(m-1,1),-1)+diag(ones(m-2,1),-2));
DD_4(1:2,1:5)=[d4;d4];DD_4(m-1:m,m-4:m)=[d4;d4];
S=-abs(q4)*(DD_4'*DD_4)/5;