-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSBP_I10.m
93 lines (65 loc) · 11.8 KB
/
SBP_I10.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
function [S,MM,BD,QQ,H,x,h] = SBP_I10(m,L)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% 10:e ordn. Pade . SBP Finita differens %%%
%%% operatorer framtagna av Ken Mattsson %%%
%%% %%%
%%% Datum: 2018 06 18 %%%
%%% %%%
%%% 7 randpunkter, bandad norm %%%
%%% %%%
%%% %%%
%%% H (Normen) %%%
%%% D1 (approx första derivatan) %%%
%%% D2 (approx andra derivatan) %%%
%%% HI*S Artificiell Dissipation %%%
%%% %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Difference (compared to Pade10)is that we optimised D1 for
% first 2 leading errors.
% For D2 we added one more Bondary point and optimised
% for first 2 leading errors.
% D1 har noggrannhet 6-10-6
% D2 har noggrannhet 4-10-4
h=L/(m-1);
x=linspace(0,L,m);
H_U=[0.2905284586765359525628226021e28 / 0.16520906281978844037605376000e29 0.314294708844475498500962467e27 / 0.1180064734427060288400384000e28 -0.196315154077194881701768309e27 / 0.660836251279153761504215040e27 0.238387881777126340635219743e27 / 0.826045314098942201880268800e27 -0.650365979101984840940006977e27 / 0.3304181256395768807521075200e28 0.743233368819318317542872677e27 / 0.8260453140989422018802688000e28 -0.17813410744959203323871303e26 / 0.786709822951373525600256000e27; 0.314294708844475498500962467e27 / 0.1180064734427060288400384000e28 0.2162217584388011375381011679e28 / 0.1376742190164903669800448000e28 -0.20651822910738607675744583e26 / 0.26223660765045784186675200e26 0.55194474455888549392089091e26 / 0.82604531409894220188026880e26 -0.175781600604149690849581787e27 / 0.550696876065961467920179200e27 0.3615283878158244084316679e25 / 0.152971354462767074422272000e27 0.318098924029259859557688541e27 / 0.8260453140989422018802688000e28; -0.196315154077194881701768309e27 / 0.660836251279153761504215040e27 -0.20651822910738607675744583e26 / 0.26223660765045784186675200e26 0.2408834556924024234358366619e28 / 0.1101393752131922935840358400e28 -0.132190815700876608783347921e27 / 0.118006473442706028840038400e27 0.78516779551706411459818891e26 / 0.122377083570213659537817600e27 -0.59527961870653558241037283e26 / 0.550696876065961467920179200e27 -0.41315042635456926196942621e26 / 0.660836251279153761504215040e27; 0.238387881777126340635219743e27 / 0.826045314098942201880268800e27 0.55194474455888549392089091e26 / 0.82604531409894220188026880e26 -0.132190815700876608783347921e27 / 0.118006473442706028840038400e27 0.145808259692653752142472099e27 / 0.68837109508245183490022400e26 -0.81111396498452467708681801e26 / 0.118006473442706028840038400e27 0.2248839138279394470409927e25 / 0.16520906281978844037605376e26 0.57640199713651591520138743e26 / 0.826045314098942201880268800e27; -0.650365979101984840940006977e27 / 0.3304181256395768807521075200e28 -0.175781600604149690849581787e27 / 0.550696876065961467920179200e27 0.78516779551706411459818891e26 / 0.122377083570213659537817600e27 -0.81111396498452467708681801e26 / 0.118006473442706028840038400e27 0.1570549637501859491337381259e28 / 0.1101393752131922935840358400e28 -0.1629918449806373679560191e25 / 0.26223660765045784186675200e26 -0.244343226922836484601858281e27 / 0.3304181256395768807521075200e28; 0.743233368819318317542872677e27 / 0.8260453140989422018802688000e28 0.3615283878158244084316679e25 / 0.152971354462767074422272000e27 -0.59527961870653558241037283e26 / 0.550696876065961467920179200e27 0.2248839138279394470409927e25 / 0.16520906281978844037605376e26 -0.1629918449806373679560191e25 / 0.26223660765045784186675200e26 0.1290633851420953729958888863e28 / 0.1376742190164903669800448000e28 0.116770637091047631171897227e27 / 0.1180064734427060288400384000e28; -0.17813410744959203323871303e26 / 0.786709822951373525600256000e27 0.318098924029259859557688541e27 / 0.8260453140989422018802688000e28 -0.41315042635456926196942621e26 / 0.660836251279153761504215040e27 0.57640199713651591520138743e26 / 0.826045314098942201880268800e27 -0.244343226922836484601858281e27 / 0.3304181256395768807521075200e28 0.116770637091047631171897227e27 / 0.1180064734427060288400384000e28 0.14619941011442378292735511573e29 / 0.16520906281978844037605376000e29;];
h4=-1/630;h3=4/315;h2=-2/45;h1=4/45;h0=8/9;
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:7,1:7)=H_U;
H(m-6:m,m-6:m)=rot90(H_U,2);
%HI=1/h*inv(H);
H=H*h;
% First derivative SBP operator, 4th order accurate at first 5 boundary points
q4=-1/280;q3=4/105;q2=-1/5;q1=4/5;
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.1197963901857751233189919e25 / 0.1259215417833753356524800e25 -0.480782067464452854859711e24 / 0.503686167133501342609920e24 0.1856654308089521347217e22 / 0.2072782580796301821440e22 -0.1466001181916337627808943e25 / 0.2518430835667506713049600e25 0.589657414614514051378571e24 / 0.2518430835667506713049600e25 -0.6238514754801649777487e22 / 0.139912824203750372947200e24; -0.1197963901857751233189919e25 / 0.1259215417833753356524800e25 0 0.57587321639578687607653e23 / 0.31091738711944527321600e23 -0.265024904459971462018211e24 / 0.157401927229219169565600e24 0.4852025514146399406691e22 / 0.3886467338993065915200e22 -0.442080975858692520107e21 / 0.728712626061199859100e21 0.71094508333989339281069e23 / 0.503686167133501342609920e24; 0.480782067464452854859711e24 / 0.503686167133501342609920e24 -0.57587321639578687607653e23 / 0.31091738711944527321600e23 0 0.164136188261950395453713e24 / 0.100737233426700268521984e24 -0.2425782872975976901711e22 / 0.1865504322716671639296e22 0.82192987107186996055e20 / 0.103639129039815091072e21 -0.565224699252225651248137e24 / 0.2518430835667506713049600e25; -0.1856654308089521347217e22 / 0.2072782580796301821440e22 0.265024904459971462018211e24 / 0.157401927229219169565600e24 -0.164136188261950395453713e24 / 0.100737233426700268521984e24 0 0.707619673683636375383459e24 / 0.503686167133501342609920e24 -0.5112295664931244065169e22 / 0.6296077089168766782624e22 0.70517279212663371878081e23 / 0.279825648407500745894400e24; 0.1466001181916337627808943e25 / 0.2518430835667506713049600e25 -0.4852025514146399406691e22 / 0.3886467338993065915200e22 0.2425782872975976901711e22 / 0.1865504322716671639296e22 -0.707619673683636375383459e24 / 0.503686167133501342609920e24 0 0.98369196616868612953181e23 / 0.93275216135833581964800e23 -0.801511293973783748447377e24 / 0.2518430835667506713049600e25; -0.589657414614514051378571e24 / 0.2518430835667506713049600e25 0.442080975858692520107e21 / 0.728712626061199859100e21 -0.82192987107186996055e20 / 0.103639129039815091072e21 0.5112295664931244065169e22 / 0.6296077089168766782624e22 -0.98369196616868612953181e23 / 0.93275216135833581964800e23 0 0.1043452766049249763128863e25 / 0.1259215417833753356524800e25; 0.6238514754801649777487e22 / 0.139912824203750372947200e24 -0.71094508333989339281069e23 / 0.503686167133501342609920e24 0.565224699252225651248137e24 / 0.2518430835667506713049600e25 -0.70517279212663371878081e23 / 0.279825648407500745894400e24 0.801511293973783748447377e24 / 0.2518430835667506713049600e25 -0.1043452766049249763128863e25 / 0.1259215417833753356524800e25 0;];
Q(1:7,1:7)=Q_U;
Q(m-6:m,m-6: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
m5=2/1575;m4=-11/1008;m3=2/63;m2=1/21;m1=-4/3;m0=4549/1800;
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 7 boundary points There was one free parameter left,
% optimised for leading order arror
M_U=[0.11648283661206993995431807e26 / 0.7915068340669306812441600e25 -0.20447863973673208563080633e26 / 0.7915068340669306812441600e25 0.9443487055341743633276417e25 / 0.3957534170334653406220800e25 -0.8971631825514136670772047e25 / 0.3957534170334653406220800e25 0.11387354646033651370545023e26 / 0.7915068340669306812441600e25 -0.4228128539923046127415177e25 / 0.7915068340669306812441600e25 0.483780379653052360771e21 / 0.5496575236575907508640e22; -0.20447863973673208563080633e26 / 0.7915068340669306812441600e25 0.26721765982109590124371e23 / 0.3435359522859942192900e22 -0.3063441525846631681284959e25 / 0.293150679284048400460800e24 0.187073864838954075595003e24 / 0.19787670851673267031104e23 -0.5413589913436301246780803e25 / 0.879452037852145201382400e24 0.2056889898002813762407e22 / 0.872472259773953572800e21 -0.634701104116935431661307e24 / 0.1583013668133861362488320e25; 0.9443487055341743633276417e25 / 0.3957534170334653406220800e25 -0.3063441525846631681284959e25 / 0.293150679284048400460800e24 0.16543432341606128487551621e26 / 0.879452037852145201382400e24 -0.2984932972470868417695355e25 / 0.158301366813386136248832e24 0.3129019036448766456733e22 / 0.261741677932186071840e21 -0.4090430838442756542383981e25 / 0.879452037852145201382400e24 0.6363995639605347909635473e25 / 0.7915068340669306812441600e25; -0.8971631825514136670772047e25 / 0.3957534170334653406220800e25 0.187073864838954075595003e24 / 0.19787670851673267031104e23 -0.2984932972470868417695355e25 / 0.158301366813386136248832e24 0.5393435737699788167189e22 / 0.245382823061424442350e21 -0.2343419024593518160655717e25 / 0.158301366813386136248832e24 0.2678880420200091099247889e25 / 0.494691771291831675777600e24 -0.722510843649482722916957e24 / 0.791506834066930681244160e24; 0.11387354646033651370545023e26 / 0.7915068340669306812441600e25 -0.5413589913436301246780803e25 / 0.879452037852145201382400e24 0.3129019036448766456733e22 / 0.261741677932186071840e21 -0.2343419024593518160655717e25 / 0.158301366813386136248832e24 0.10477021842714816355940011e26 / 0.879452037852145201382400e24 -0.495846774760156412050523e24 / 0.97716893094682800153600e23 0.2789907087983859616938271e25 / 0.3957534170334653406220800e25; -0.4228128539923046127415177e25 / 0.7915068340669306812441600e25 0.2056889898002813762407e22 / 0.872472259773953572800e21 -0.4090430838442756542383981e25 / 0.879452037852145201382400e24 0.2678880420200091099247889e25 / 0.494691771291831675777600e24 -0.495846774760156412050523e24 / 0.97716893094682800153600e23 0.110156959340418440231647e24 / 0.27482876182879537543200e23 -0.12593657480195191556140951e26 / 0.7915068340669306812441600e25; 0.483780379653052360771e21 / 0.5496575236575907508640e22 -0.634701104116935431661307e24 / 0.1583013668133861362488320e25 0.6363995639605347909635473e25 / 0.7915068340669306812441600e25 -0.722510843649482722916957e24 / 0.791506834066930681244160e24 0.2789907087983859616938271e25 / 0.3957534170334653406220800e25 -0.12593657480195191556140951e26 / 0.7915068340669306812441600e25 0.4070677235094950940663229e25 / 0.1583013668133861362488320e25;];
M(1:7,1:7)=M_U;
M(m-6:m,m-6:m)=rot90(M_U,2);
M=M/h;
d1_U=[-0.49e2 / 0.20e2 6 -0.15e2 / 0.2e1 0.20e2 / 0.3e1 -0.15e2 / 0.4e1 0.6e1 / 0.5e1 -0.1e1 / 0.6e1;]/h;
d1_l=zeros(1,m);
d1_l(1:7)=d1_U;
d1_r=zeros(1,m);
d1_r(m-6: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 5-11-5
d6=0*[1 -6 15 -20 15 -6 1];
DD_6=(diag(ones(m-3,1),3)-6*diag(ones(m-2,1),2)+15*diag(ones(m-1,1),1)-20*diag(ones(m,1),0)+15*diag(ones(m-1,1),-1)-6*diag(ones(m-2,1),-2)+diag(ones(m-3,1),-3));
DD_6(1:3,1:7)=[d6;d6;d6];DD_6(m-2:m,m-6:m)=[d6;d6;d6];
DD=DD_6'*DD_6;
aa=max(max(DD));
S=-1/aa*DD;