Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_compatible_2.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_2(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);H(1,1)=1/2;H(m,m)=1/2;
30
31
32 H=H*h;
33 HI=inv(H);
34
35
36 % First derivative SBP operator, 1st order accurate at first 6 boundary points
37
38 q1=1/2;
39 Q=q1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1));
40
41 %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));
42
43
44 e_1=zeros(m,1);e_1(1)=1;
45 e_m=zeros(m,1);e_m(m)=1;
46
47
48 D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ;
49
50 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51
52
53
54 % Second derivative, 1st order accurate at first 6 boundary points
55 m1=-1;m0=2;
56 M=m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0);M(1,1)=1;M(m,m)=1;
57 M=M/h;
58
59 S_U=[-1 1]/h;
60 S_1=zeros(1,m);
61 S_1(1:2)=S_U;
62 S_m=zeros(1,m);
63
64 S_m(m-1:m)=fliplr(-S_U);
65
66 D2=HI*(-M-e_1*S_1+e_m*S_m);
67
68
69 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
70
71
72
73 % Third derivative, 1st order accurate at first 6 boundary points
74
75 q2=1/2;q1=-1;
76 Q3=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));
77
78 %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));
79
80
81 Q3_U = [0 -0.2e1 / 0.5e1 0.3e1 / 0.10e2 0.1e1 / 0.10e2; 0.2e1 / 0.5e1 0 -0.7e1 / 0.10e2 0.3e1 / 0.10e2; -0.3e1 / 0.10e2 0.7e1 / 0.10e2 0 -0.9e1 / 0.10e2; -0.1e1 / 0.10e2 -0.3e1 / 0.10e2 0.9e1 / 0.10e2 0;];
82 Q3(1:4,1:4)=Q3_U;
83 Q3(m-3:m,m-3:m)=flipud( fliplr( -Q3_U ) );
84 Q3=Q3/h^2;
85
86
87
88 S2_U=[1 -2 1;]/h^2;
89 S2_1=zeros(1,m);
90 S2_1(1:3)=S2_U;
91 S2_m=zeros(1,m);
92 S2_m(m-2:m)=fliplr(S2_U);
93
94
95
96 D3=HI*(Q3 - e_1*S2_1 + e_m*S2_m +1/2*S_1'*S_1 -1/2*S_m'*S_m ) ;
97
98 % Fourth derivative, 0th order accurate at first 6 boundary points (still
99 % yield 4th order convergence if stable: for example u_tt=-u_xxxx
100
101 m2=1;m1=-4;m0=6;
102 M4=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);
103
104 %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));
105
106 M4_U=[0.4e1 / 0.5e1 -0.7e1 / 0.5e1 0.2e1 / 0.5e1 0.1e1 / 0.5e1; -0.7e1 / 0.5e1 0.16e2 / 0.5e1 -0.11e2 / 0.5e1 0.2e1 / 0.5e1; 0.2e1 / 0.5e1 -0.11e2 / 0.5e1 0.21e2 / 0.5e1 -0.17e2 / 0.5e1; 0.1e1 / 0.5e1 0.2e1 / 0.5e1 -0.17e2 / 0.5e1 0.29e2 / 0.5e1;];
107
108 M4(1:4,1:4)=M4_U;
109
110 M4(m-3:m,m-3:m)=flipud( fliplr( M4_U ) );
111 M4=M4/h^3;
112
113 S3_U=[-1 3 -3 1;]/h^3;
114 S3_1=zeros(1,m);
115 S3_1(1:4)=S3_U;
116 S3_m=zeros(1,m);
117 S3_m(m-3:m)=fliplr(-S3_U);
118
119 D4=HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m);
120
121
122
123 S_1 = S_1';
124 S_m = S_m';
125 S2_1 = S2_1';
126 S2_m = S2_m';
127 S3_1 = S3_1';
128 S3_m = S3_m';
129
130
131
132
133 % L=h*(m-1);
134
135 % x1=linspace(0,L,m)';
136 % x2=x1.^2/fac(2);
137 % x3=x1.^3/fac(3);
138 % x4=x1.^4/fac(4);
139 % x5=x1.^5/fac(5);
140
141 % x0=x1.^0/fac(1);
142
143
144 end