Mercurial > repos > public > sbplib
annotate +sbp/+implementations/d4_variable_6_3.m @ 318:99005a80b4c2 feature/beams
Cleaned up d4_variable_4_min_boundary. Removed incorrect D2s from a bunch of files.
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Mon, 26 Sep 2016 08:44:17 +0200 |
parents | 203afa156f59 |
children | def409c10800 |
rev | line source |
---|---|
312 | 1 function [H, HI, D1, D2, D4, e_l, e_r, M4, d2_l, d2_r, d3_l, d3_r, d1_l, d1_r] = d4_variable_6_3(m,h) |
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
3 %%% 6:te ordn. SBP Finita differens %%% | |
4 %%% operatorer med diagonal norm %%% | |
5 %%% Extension to variable koeff %%% | |
6 %%% %%% | |
7 %%% H (Normen) %%% | |
8 %%% D1=H^(-1)Q (approx f?rsta derivatan) %%% | |
9 %%% D2 (approx andra derivatan) %%% | |
10 %%% D2=HI*(R+C*D*S %%% | |
11 %%% %%% | |
12 %%% R=-D1'*H*C*D1-RR %%% | |
13 %%% %%% | |
14 %%% RR ?r dissipation) %%% | |
15 %%% Dissipationen uppbyggd av D4: %%% | |
16 %%% DI=D4*B*H*D4 %%% | |
17 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
18 |
312 | 19 % H?r med 7 RP ist?llet f?r 8 f?r D4 operatorn, dock samma randderivator |
20 % Denna ?r noggrannare, och har 2a ordningens randdslutning och b?r ge 6te | |
21 % ordningens konvergens. Hade 2 fria parametrar att optimera | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
22 |
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
23 |
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
24 |
312 | 25 H=diag(ones(m,1),0); |
26 H(1:7,1:7)=[ | |
27 0.414837907e9/0.1191965760e10 0 0 0 0 0 0; | |
28 0 0.475278367e9/0.397321920e9 0 0 0 0 0; | |
29 0 0 0.13872751e8/0.12416310e8 0 0 0 0; | |
30 0 0 0 0.346739027e9/0.595982880e9 0 0 0; | |
31 0 0 0 0 0.560227469e9/0.397321920e9 0 0; | |
32 0 0 0 0 0 0.322971631e9/0.397321920e9 0; | |
33 0 0 0 0 0 0 0.616122491e9/0.595982880e9; | |
34 ]; | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
35 |
312 | 36 H(m-6:m,m-6:m) = fliplr(flipud(H(1:7,1:7))); |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
37 |
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
38 |
312 | 39 e_1=zeros(m,1); |
40 e_1(1)=1; | |
41 e_m=zeros(m,1); | |
42 e_m(m)=1; | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
43 |
312 | 44 S_U=[-0.137e3/0.60e2 5 -5 0.10e2/0.3e1 -0.5e1/0.4e1 0.1e1/0.5e1;]/h; |
45 S_1=zeros(1,m); | |
46 S_1(1:6)=S_U; | |
47 S_m=zeros(1,m); | |
48 S_m(m-5:m)=fliplr(-S_U); | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
49 |
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
50 |
316
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
51 S2_U = [0.15e2/0.4e1 -0.77e2/0.6e1 0.107e3/0.6e1 -13 0.61e2/0.12e2 -0.5e1/0.6e1;]/h^2; |
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
52 S2_1 = zeros(1,m); |
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
53 S2_1(1:6) = S2_U; |
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
54 S2_m = zeros(1,m); |
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
55 S2_m(m-5:m) = fliplr(S2_U); |
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
56 |
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
57 |
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
58 S3_U = [-0.17e2/0.4e1 0.71e2/0.4e1 -0.59e2/0.2e1 0.49e2/0.2e1 -0.41e2/0.4e1 0.7e1/0.4e1;]/h^3; |
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
59 S3_1 = zeros(1,m); |
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
60 S3_1(1:6) = S3_U; |
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
61 S3_m = zeros(1,m); |
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
62 S3_m(m-5:m) = fliplr(-S3_U); |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
63 |
312 | 64 %DS=zeros(m,m); |
65 %DS(1,1:5)=-[-25/12, 4, -3, 4/3, -1/4]; | |
66 %DS(m,m-4:m)=fliplr(-[-25/12, 4, -3, 4/3, -1/4]); | |
67 %DS=diag(c)*DS/h; | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
68 |
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
69 |
312 | 70 H=h*H; |
71 HI=inv(H); | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
72 |
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
73 |
312 | 74 % Fourth derivative, 1th order accurate at first 8 boundary points (still |
75 % yield 5th order convergence if stable: for example u_tt=-u_xxxx | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
76 |
312 | 77 m4 = 7/240; |
78 m3 = -2/5; | |
79 m2 = 169/60; | |
80 m1 = -122/15; | |
81 m0 = 91/8; | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
82 |
312 | 83 M4 = 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); |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
84 |
312 | 85 %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)); |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
86 |
312 | 87 M4_U = [ |
88 0.1399708478939e13/0.263487168000e12 -0.13482796013041e14/0.834376032000e12 0.344344095859e12/0.17565811200e11 -0.3166261424681e13/0.250312809600e12 0.1508605165681e13/0.333750412800e12 -0.486270829441e12/0.834376032000e12 -0.221976356359e12/0.5006256192000e13; | |
89 -0.13482796013041e14/0.834376032000e12 0.7260475818391e13/0.139062672000e12 -0.27224036353e11/0.406022400e9 0.1847477458951e13/0.41718801600e11 -0.848984558161e12/0.55625068800e11 0.247494925991e12/0.139062672000e12 0.165585445559e12/0.834376032000e12; | |
90 0.344344095859e12/0.17565811200e11 -0.27224036353e11/0.406022400e9 0.2044938640393e13/0.22250027520e11 -0.1071086785417e13/0.16687520640e11 0.502199537033e12/0.22250027520e11 -0.143589154441e12/0.55625068800e11 -0.88181965559e11/0.333750412800e12; | |
91 -0.3166261424681e13/0.250312809600e12 0.1847477458951e13/0.41718801600e11 -0.1071086785417e13/0.16687520640e11 0.628860435593e12/0.12515640480e11 -0.73736245829e11/0.3337504128e10 0.195760572271e12/0.41718801600e11 -0.81156046361e11/0.250312809600e12; | |
92 0.1508605165681e13/0.333750412800e12 -0.848984558161e12/0.55625068800e11 0.502199537033e12/0.22250027520e11 -0.73736245829e11/0.3337504128e10 0.76725285869e11/0.4450005504e10 -0.3912429433e10/0.406022400e9 0.53227370659e11/0.17565811200e11; | |
93 -0.486270829441e12/0.834376032000e12 0.247494925991e12/0.139062672000e12 -0.143589154441e12/0.55625068800e11 0.195760572271e12/0.41718801600e11 -0.3912429433e10/0.406022400e9 0.1699707221791e13/0.139062672000e12 -0.6959018412841e13/0.834376032000e12; | |
94 -0.221976356359e12/0.5006256192000e13 0.165585445559e12/0.834376032000e12 -0.88181965559e11/0.333750412800e12 -0.81156046361e11/0.250312809600e12 0.53227370659e11/0.17565811200e11 -0.6959018412841e13/0.834376032000e12 0.3012195053939e13/0.263487168000e12; | |
95 ]; | |
96 | |
97 M4(1:7,1:7) = M4_U; | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
98 |
312 | 99 M4(m-6:m,m-6:m) = flipud( fliplr( M4_U ) ); |
100 M4 = M4/h^3; | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
101 |
312 | 102 D4 = HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m); |
103 end |