-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSBP_SL6.m
115 lines (83 loc) · 9.79 KB
/
SBP_SL6.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_SL6(m,L)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 6:te Spectral SBP Finita differens %%%
%%% operatorer framtagna av Ken Mattsson %%%
%%% %%%
%%% Datum: 2019 07 10 %%%
%%% %%%
%%% 5 randpunkter, %%%
%%% %%%
%%% Opt. för w=1.4, 1.8, 2.1, 2.4, 2.6 %%%
%%% %%%
%%% Neptadiagonal norm %%%
%%% Heptadiagonal Q %%%
%%% %%%
%%% Noggrannhet 4-6-4 %%%
%%% %%%
%%% %%%
%%% %%%
%%% H (Normen) %%%
%%% D1 (approx första derivatan) %%%
%%% D2 (approx andra derivatan) %%%
%%% S Artificiell Dissipation %%%
%%% %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h=L/(m-1);
x=linspace(0,L,m)';
% D1 har noggrannhet 3-4-3
H_U=[0.79878580442295032441601874486737144415280043647933755752713e59 / 0.402067784072662543662816727623787236373357948281439208976000e60 0.1042403614517472619442893571774261837350269541915086316151e58 / 0.4730209224384265219562549736750438074980681744487520105600e58 -0.64858110001848049165080897384691266778758048947270846454853e59 / 0.402067784072662543662816727623787236373357948281439208976000e60 0.5345370115598851919154692205566064498836549400913491079513e58 / 0.44674198230295838184757414180420804041484216475715467664000e59 -0.1e1 / 0.40e2; 0.1042403614517472619442893571774261837350269541915086316151e58 / 0.4730209224384265219562549736750438074980681744487520105600e58 0.532927202702082527103745280796986059837192710581206481381713e60 / 0.402067784072662543662816727623787236373357948281439208976000e60 -0.5186430119890140107466611920844281189498631850909806434079e58 / 0.16082711362906501746512669104951489454934317931257568359040e59 0.1799568730157969765512867074872430793674221397852795443387e58 / 0.402067784072662543662816727623787236373357948281439208976000e60 -0.1e1 / 0.125e3; -0.64858110001848049165080897384691266778758048947270846454853e59 / 0.402067784072662543662816727623787236373357948281439208976000e60 -0.5186430119890140107466611920844281189498631850909806434079e58 / 0.16082711362906501746512669104951489454934317931257568359040e59 0.61256111037224830190998678820219373216896503501682222451017e59 / 0.44674198230295838184757414180420804041484216475715467664000e59 0.1034346157536012084946870188918571623876925590398091330523e58 / 0.16082711362906501746512669104951489454934317931257568359040e59 -0.199658573871448799093067904895417890095163280294709171427e57 / 0.6282309126135352244731511369121675568333717941897487640250e58; 0.5345370115598851919154692205566064498836549400913491079513e58 / 0.44674198230295838184757414180420804041484216475715467664000e59 0.1799568730157969765512867074872430793674221397852795443387e58 / 0.402067784072662543662816727623787236373357948281439208976000e60 0.1034346157536012084946870188918571623876925590398091330523e58 / 0.16082711362906501746512669104951489454934317931257568359040e59 0.245321822725843042794132332652334614326573095845596171448313e60 / 0.402067784072662543662816727623787236373357948281439208976000e60 0.345547106745479613454633339937488453986601951173818465773e57 / 0.1256461825227070448946302273824335113666743588379497528050e58; -0.1e1 / 0.40e2 -0.1e1 / 0.125e3 -0.199658573871448799093067904895417890095163280294709171427e57 / 0.6282309126135352244731511369121675568333717941897487640250e58 0.345547106745479613454633339937488453986601951173818465773e57 / 0.1256461825227070448946302273824335113666743588379497528050e58 0.14888551080199451729258667644794476158622993569253663573893e59 / 0.25129236504541408978926045476486702273334871767589950561000e59;];
h0=4203267613564094932432577824954/7049220443079284250976145948443;
h1=22618790744689935699264926210401/84590645316951411011713751381316;
h2=-2209778222820418388602425303685/42295322658475705505856875690658;
h3=-1581945765/75409415044;
h4=228992488/33235651987;
h5=27214243/33751459947;
H=h5*(diag(ones(m-5,1),5)+diag(ones(m-5,1),-5))+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:5,1:5)=H_U;
H(m-4:m,m-4:m)=rot90(H_U,2);
H=H*h;
q1=9607266784889201296177/19560081711822931675052;
q2=8866705546306148289391/97800408559114658375260;
q3=-19659090145677941034997/293401225677343975125780;
q4=127051314/37983174851;
q5=389910724/128741750713;
Q=q5*(diag(ones(m-5,1),5) - diag(ones(m-5,1),-5))+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.15239128396755400249969206956165016990021961372563121379e56 / 0.18614249262623265910315589241842001683951756864881444860e56 -0.8954259994747811097147818923183180013611755023851765877813e58 / 0.16752824336360939319284030317657801515556581178393300374000e59 0.4402359768928196260518664819276624736003210139681706890413e58 / 0.16752824336360939319284030317657801515556581178393300374000e59 -0.1e1 / 0.20e2; -0.15239128396755400249969206956165016990021961372563121379e56 / 0.18614249262623265910315589241842001683951756864881444860e56 0 0.4370906932242860734619186336277091246058287649946015161813e58 / 0.3350564867272187863856806063531560303111316235678660074800e58 -0.3085067519172593588524174237810672596749725053718816631551e58 / 0.5584274778786979773094676772552600505185527059464433458000e58 0.63069266460854237084961553931969421315853888745750096168e56 / 0.1047051521022558707455251894853612594722286323649581273375e58; 0.8954259994747811097147818923183180013611755023851765877813e58 / 0.16752824336360939319284030317657801515556581178393300374000e59 -0.4370906932242860734619186336277091246058287649946015161813e58 / 0.3350564867272187863856806063531560303111316235678660074800e58 0 0.342493552108155207944860764638043186526514517135856219021e57 / 0.418820608409023482982100757941445037888914529459832509350e57 0.27033199915418931111584781689776447769064409521277839689e56 / 0.2094103042045117414910503789707225189444572647299162546750e58; -0.4402359768928196260518664819276624736003210139681706890413e58 / 0.16752824336360939319284030317657801515556581178393300374000e59 0.3085067519172593588524174237810672596749725053718816631551e58 / 0.5584274778786979773094676772552600505185527059464433458000e58 -0.342493552108155207944860764638043186526514517135856219021e57 / 0.418820608409023482982100757941445037888914529459832509350e57 0 0.4635444786847164833869474111550667174723773195694372071e55 / 0.9307124631311632955157794620921000841975878432440722430e55; 0.1e1 / 0.20e2 -0.63069266460854237084961553931969421315853888745750096168e56 / 0.1047051521022558707455251894853612594722286323649581273375e58 -0.27033199915418931111584781689776447769064409521277839689e56 / 0.2094103042045117414910503789707225189444572647299162546750e58 -0.4635444786847164833869474111550667174723773195694372071e55 / 0.9307124631311632955157794620921000841975878432440722430e55 0;];
Q(1:5,1:5)=Q_U;
Q(m-4:m,m-4: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, 3rd order accurate at first 6 boundary points
% This one optimized dispersion with 2 adional free coefficients (m5 and
% m5), chosen to be exactat w=2 and w=2.6
m0 = 0.85086800469647609167;
m1 = -0.078209094071890664434;
m2 = -0.42716031111161651324;
m3 = 0.059459556975949098755;
m4 = 0.028723540890945572833;
m5 = -0.0080958811471095540142;
m6 = -0.00015181388451598573244;
M=m6*(diag(ones(m-6,1),6)+diag(ones(m-6,1),-6))+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 error
M_U=[0.13287327019037196101e1 -0.18900389300917004448e1 0.85399456087430111179e0 -0.36211353131161676938e0 0.78583776964017946454e-1 -0.90067644542054684088e-2; -0.18900389300917004448e1 0.45758916623385576471e1 -0.36207685103085648816e1 0.10883712991791768949e1 -0.17620857524452903462e0 0.31000749158685358820e-1; 0.85399456087430111179e0 -0.36207685103085648816e1 0.47392083449318231267e1 -0.19431139403470281309e1 -0.10621958029548064316e0 0.56423279285629384107e-1; -0.36211353131161676938e0 0.10883712991791768949e1 -0.19431139403470281309e1 0.18166321481728326986e1 -0.25482827568475709780e0 -0.42488310284387672726e0; 0.78583776964017946454e-1 -0.17620857524452903462e0 -0.10621958029548064316e0 -0.25482827568475709780e0 0.88501753991608278936e0 -0.79119977378986578829e-1; -0.90067644542054684088e-2 0.31000749158685358820e-1 0.56423279285629384107e-1 -0.42488310284387672726e0 -0.79119977378986578829e-1 0.85101981858099207740e0;];
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];
S=-abs(q5)*DD_5'*DD_5;