Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_compatible_4.m @ 261:6009f2712d13 operator_remake

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Moved and renamned all implementations.
author Martin Almquist <martin.almquist@it.uu.se>
date Thu, 08 Sep 2016 15:35:45 +0200
parents
children bfa130b7abf6
comparison
equal deleted inserted replaced
260:b4116ce49ac4 261:6009f2712d13
1 function [H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m,...
2 S3_1, S3_m, S_1, S_m] = d4_compatible_4(m,h)
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 %%% 4:de ordn. SBP Finita differens %%%
5 %%% operatorer framtagna av Ken Mattsson %%%
6 %%% %%%
7 %%% 6 randpunkter, diagonal norm %%%
8 %%% %%%
9 %%% Datum: 2013-11-11 %%%
10 %%% %%%
11 %%% %%%
12 %%% H (Normen) %%%
13 %%% D1 (approx f?rsta derivatan) %%%
14 %%% D2 (approx andra derivatan) %%%
15 %%% D3 (approx tredje derivatan) %%%
16 %%% D2 (approx fj?rde derivatan) %%%
17 %%% %%%
18 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
19
20 % M?ste ange antal punkter (m) och stegl?ngd (h)
21 % Notera att dessa opetratorer ?r framtagna f?r att anv?ndas n?r
22 % vi har 3de och 4de derivator i v?r PDE
23 % I annat fall anv?nd de "traditionella" som har noggrannare
24 % randsplutningar f?r D1 och D2
25
26 % Vi b?rjar med normen. Notera att alla SBP operatorer delar samma norm,
27 % vilket ?r n?dv?ndigt f?r stabilitet
28
29 H=diag(ones(m,1),0);
30 H_U=[0.3e1 / 0.11e2 0 0 0 0 0; 0 0.2125516311e10 / 0.1311004640e10 0 0 0 0; 0 0 0.278735189e9 / 0.1966506960e10 0 0 0; 0 0 0 0.285925927e9 / 0.163875580e9 0 0; 0 0 0 0 0.1284335339e10 / 0.1966506960e10 0; 0 0 0 0 0 0.4194024163e10 / 0.3933013920e10;];
31 H(1:6,1:6)=H_U;
32 H(m-5:m,m-5:m)=fliplr(flipud(H_U));
33 H=H*h;
34 HI=inv(H);
35
36
37 % First derivative SBP operator, 1st order accurate at first 6 boundary points
38
39 q2=-1/12;q1=8/12;
40 Q=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));
41
42 %Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2));
43
44 Q_U = [0 0.9e1 / 0.11e2 -0.9e1 / 0.22e2 0.1e1 / 0.11e2 0 0; -0.9e1 / 0.11e2 0 0.2595224893e10 / 0.2622009280e10 -0.151435707e9 / 0.327751160e9 0.1112665611e10 / 0.2622009280e10 -0.1290899e7 / 0.9639740e7; 0.9e1 / 0.22e2 -0.2595224893e10 / 0.2622009280e10 0 0.1468436423e10 / 0.983253480e9 -0.1194603401e10 / 0.983253480e9 0.72033031e8 / 0.238364480e9; -0.1e1 / 0.11e2 0.151435707e9 / 0.327751160e9 -0.1468436423e10 / 0.983253480e9 0 0.439819541e9 / 0.327751160e9 -0.215942641e9 / 0.983253480e9; 0 -0.1112665611e10 / 0.2622009280e10 0.1194603401e10 / 0.983253480e9 -0.439819541e9 / 0.327751160e9 0 0.1664113643e10 / 0.2622009280e10; 0 0.1290899e7 / 0.9639740e7 -0.72033031e8 / 0.238364480e9 0.215942641e9 / 0.983253480e9 -0.1664113643e10 / 0.2622009280e10 0;];
45 Q(1:6,1:6)=Q_U;
46 Q(m-5:m,m-5:m)=flipud( fliplr( -Q_U ) );
47
48 e_1=zeros(m,1);e_1(1)=1;
49 e_m=zeros(m,1);e_m(m)=1;
50
51
52 D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ;
53
54 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
55
56
57
58 % % Second derivative, 1st order accurate at first 6 boundary points
59 % m2=1/12;m1=-16/12;m0=30/12;
60 % M=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);
61 % %M=(1/12*diag(ones(m-2,1),2)-16/12*diag(ones(m-1,1),1)-16/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2)+30/12*diag(ones(m,1),0));
62 % M_U=[0.2386127e7 / 0.2177280e7 -0.515449e6 / 0.453600e6 -0.10781e5 / 0.777600e6 0.61567e5 / 0.1360800e7 0.6817e4 / 0.403200e6 -0.1069e4 / 0.136080e6; -0.515449e6 / 0.453600e6 0.4756039e7 / 0.2177280e7 -0.1270009e7 / 0.1360800e7 -0.3751e4 / 0.28800e5 0.3067e4 / 0.680400e6 0.119459e6 / 0.10886400e8; -0.10781e5 / 0.777600e6 -0.1270009e7 / 0.1360800e7 0.111623e6 / 0.60480e5 -0.555587e6 / 0.680400e6 -0.551339e6 / 0.5443200e7 0.8789e4 / 0.453600e6; 0.61567e5 / 0.1360800e7 -0.3751e4 / 0.28800e5 -0.555587e6 / 0.680400e6 0.1025327e7 / 0.544320e6 -0.464003e6 / 0.453600e6 0.222133e6 / 0.5443200e7; 0.6817e4 / 0.403200e6 0.3067e4 / 0.680400e6 -0.551339e6 / 0.5443200e7 -0.464003e6 / 0.453600e6 0.5074159e7 / 0.2177280e7 -0.1784047e7 / 0.1360800e7; -0.1069e4 / 0.136080e6 0.119459e6 / 0.10886400e8 0.8789e4 / 0.453600e6 0.222133e6 / 0.5443200e7 -0.1784047e7 / 0.1360800e7 0.1812749e7 / 0.725760e6;];
63 %
64 % M(1:6,1:6)=M_U;
65 %
66 % M(m-5:m,m-5:m)=flipud( fliplr( M_U ) );
67 % M=M/h;
68 %
69 S_U=[-0.11e2 / 0.6e1 3 -0.3e1 / 0.2e1 0.1e1 / 0.3e1;]/h;
70 S_1=zeros(1,m);
71 S_1(1:4)=S_U;
72 S_m=zeros(1,m);
73
74 S_m(m-3:m)=fliplr(-S_U);
75
76 % D2=HI*(-M-e_1*S_1+e_m*S_m);
77
78
79 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
80
81
82
83 % Third derivative, 1st order accurate at first 6 boundary points
84
85 % q3=-1/8;q2=1;q1=-13/8;
86 % Q3=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));
87 %
88 % %QQ3=(-1/8*diag(ones(m-3,1),3) + 1*diag(ones(m-2,1),2) - 13/8*diag(ones(m-1,1),1) +13/8*diag(ones(m-1,1),-1) -1*diag(ones(m-2,1),-2) + 1/8*diag(ones(m-3,1),-3));
89 %
90 %
91 % Q3_U = [0 -0.88471e5 / 0.67200e5 0.58139e5 / 0.33600e5 -0.1167e4 / 0.2800e4 -0.89e2 / 0.11200e5 0.7e1 / 0.640e3; 0.88471e5 / 0.67200e5 0 -0.43723e5 / 0.16800e5 0.46783e5 / 0.33600e5 -0.191e3 / 0.3200e4 -0.1567e4 / 0.33600e5; -0.58139e5 / 0.33600e5 0.43723e5 / 0.16800e5 0 -0.4049e4 / 0.2400e4 0.29083e5 / 0.33600e5 -0.71e2 / 0.1400e4; 0.1167e4 / 0.2800e4 -0.46783e5 / 0.33600e5 0.4049e4 / 0.2400e4 0 -0.8591e4 / 0.5600e4 0.10613e5 / 0.11200e5; 0.89e2 / 0.11200e5 0.191e3 / 0.3200e4 -0.29083e5 / 0.33600e5 0.8591e4 / 0.5600e4 0 -0.108271e6 / 0.67200e5; -0.7e1 / 0.640e3 0.1567e4 / 0.33600e5 0.71e2 / 0.1400e4 -0.10613e5 / 0.11200e5 0.108271e6 / 0.67200e5 0;];
92 %
93 % Q3(1:6,1:6)=Q3_U;
94 % Q3(m-5:m,m-5:m)=flipud( fliplr( -Q3_U ) );
95 % Q3=Q3/h^2;
96
97
98
99 S2_U=[2 -5 4 -1;]/h^2;
100 S2_1=zeros(1,m);
101 S2_1(1:4)=S2_U;
102 S2_m=zeros(1,m);
103 S2_m(m-3:m)=fliplr(S2_U);
104
105
106
107 %D3=HI*(Q3 - e_1*S2_1 + e_m*S2_m +1/2*S_1'*S_1 -1/2*S_m'*S_m ) ;
108
109 % Fourth derivative, 0th order accurate at first 6 boundary points (still
110 % yield 4th order convergence if stable: for example u_tt=-u_xxxx
111
112 m3=-1/6;m2=2;m1=-13/2;m0=28/3;
113 M4=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);
114
115 %M4=(-1/6*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3) ) + 2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2)) -13/2*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1)) + 28/3*diag(ones(m,1),0));
116
117 M4_U=[0.227176919517319e15 / 0.94899692875680e14 -0.15262605263734e14 / 0.2965615402365e13 0.20205404771243e14 / 0.6778549491120e13 -0.3998303664097e13 / 0.23724923218920e14 0.1088305091927e13 / 0.94899692875680e14 -0.1686077077693e13 / 0.23724923218920e14; -0.15262605263734e14 / 0.2965615402365e13 0.280494781164181e15 / 0.23724923218920e14 -0.46417445546261e14 / 0.5931230804730e13 0.1705307929429e13 / 0.1694637372780e13 -0.553547394061e12 / 0.5931230804730e13 0.5615721694973e13 / 0.23724923218920e14; 0.20205404771243e14 / 0.6778549491120e13 -0.46417445546261e14 / 0.5931230804730e13 0.4135802350237e13 / 0.551742400440e12 -0.4140981465247e13 / 0.1078405600860e13 0.75538453067437e14 / 0.47449846437840e14 -0.4778134936391e13 / 0.11862461609460e14; -0.3998303664097e13 / 0.23724923218920e14 0.1705307929429e13 / 0.1694637372780e13 -0.4140981465247e13 / 0.1078405600860e13 0.20760974175677e14 / 0.2965615402365e13 -0.138330689701889e15 / 0.23724923218920e14 0.23711317526909e14 / 0.11862461609460e14; 0.1088305091927e13 / 0.94899692875680e14 -0.553547394061e12 / 0.5931230804730e13 0.75538453067437e14 / 0.47449846437840e14 -0.138330689701889e15 / 0.23724923218920e14 0.120223780251937e15 / 0.13557098982240e14 -0.151383731537477e15 / 0.23724923218920e14; -0.1686077077693e13 / 0.23724923218920e14 0.5615721694973e13 / 0.23724923218920e14 -0.4778134936391e13 / 0.11862461609460e14 0.23711317526909e14 / 0.11862461609460e14 -0.151383731537477e15 / 0.23724923218920e14 0.220304030094121e15 / 0.23724923218920e14;];
118
119 M4(1:6,1:6)=M4_U;
120
121 M4(m-5:m,m-5:m)=flipud( fliplr( M4_U ) );
122 M4=M4/h^3;
123
124 S3_U=[-1 3 -3 1;]/h^3;
125 S3_1=zeros(1,m);
126 S3_1(1:4)=S3_U;
127 S3_m=zeros(1,m);
128 S3_m(m-3:m)=fliplr(-S3_U);
129
130 D4=HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m);
131
132 S_1 = S_1';
133 S_m = S_m';
134 S2_1 = S2_1';
135 S2_m = S2_m';
136 S3_1 = S3_1';
137 S3_m = S3_m';
138
139 % L=h*(m-1);
140 %
141 % x1=linspace(0,L,m)';
142 % x2=x1.^2/fac(2);
143 % x3=x1.^3/fac(3);
144 % x4=x1.^4/fac(4);
145 % x5=x1.^5/fac(5);
146 %
147 % x0=x1.^0/fac(1);
148
149 end