Mercurial > repos > public > sbplib
annotate +sbp/+implementations/d4_variable_6_min_boundary_points.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 | fc91e3a2d44a |
| 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_min_boundary_points(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 6 RP ist?llet f?r 8 f?r D4 operatorn, dock samma randderivator |
|
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
20 |
| 312 | 21 H = diag(ones(m,1),0); |
| 22 H(1:6,1:6) = diag([13649/43200,12013/8640,2711/4320,5359/4320,7877/8640, 43801/43200]); | |
| 23 H(m-5:m,m-5:m) = fliplr(flipud(diag([13649/43200,12013/8640, 2711/4320,5359/4320,7877/8640,43801/43200]))); | |
|
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
24 |
|
316
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
25 e_1 = zeros(m,1);e_1(1)=1; |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
26 e_m = zeros(m,1);e_m(m)=1; |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
27 |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
28 S_U = [-25/12, 4, -3, 4/3, -1/4]/h; |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
29 S_1 = zeros(1,m); |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
30 S_1(1:5) = S_U; |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
31 S_m = zeros(1,m); |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
32 S_m(m-4:m) = fliplr(-S_U); |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
33 |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
34 S2_U = [0.35e2/0.12e2 -0.26e2/0.3e1 0.19e2/0.2e1 -0.14e2/0.3e1 0.11e2/0.12e2;]/h^2; |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
35 S2_1 = zeros(1,m); |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
36 S2_1(1:5) = S2_U; |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
37 S2_m = zeros(1,m); |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
38 S2_m(m-4:m) = fliplr(S2_U); |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
39 |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
40 S3_U = [-0.5e1/0.2e1 9 -12 7 -0.3e1/0.2e1;]/h^3; |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
41 S3_1 = zeros(1,m); |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
42 S3_1(1:5) = S3_U; |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
43 S3_m = zeros(1,m); |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
44 S3_m(m-4:m) = fliplr(-S3_U); |
|
203afa156f59
Collected boundary operators.
Jonatan Werpers <jonatan@werpers.com>
parents:
312
diff
changeset
|
45 |
|
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
46 |
| 312 | 47 x1=0.70127127127127; |
|
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
48 |
|
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
49 |
| 312 | 50 D1=(1/60*diag(ones(m-3,1),3)-9/60*diag(ones(m-2,1),2)+45/60*diag(ones(m-1,1),1)-45/60*diag(ones(m-1,1),-1)+9/60*diag(ones(m-2,1),-2)-1/60*diag(ones(m-3,1),-3)); |
|
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
51 |
| 312 | 52 D1(1:6,1:9) = [ |
| 53 -21600/13649, 43200/13649*x1-7624/40947, -172800/13649*x1+ 715489/81894, 259200/13649*x1-187917/13649, -172800/13649* x1+735635/81894, 43200/13649*x1-89387/40947, 0, 0, 0; | |
| 54 -8640/12013*x1+7624/180195, 0, 86400/12013*x1-57139/12013, -172800/12013*x1+745733/72078, 129600/12013*x1-91715/12013, -34560/12013*x1+240569/120130, 0, 0, 0; | |
| 55 17280/2711*x1-715489/162660, -43200/2711*x1+57139/5422, 0, 86400/2711*x1-176839/8133, -86400/2711*x1+242111/10844, 25920/2711*x1-182261/27110, 0, 0, 0; | |
| 56 -25920/5359*x1+187917/53590, 86400/5359*x1-745733/64308, -86400/5359*x1+176839/16077, 0, 43200/5359*x1-165041/32154, -17280/5359*x1+710473/321540, 72/5359, 0, 0; | |
| 57 34560/7877*x1-147127/47262, -129600/7877*x1+91715/7877, 172800/7877*x1-242111/15754, -86400/7877*x1+165041/23631, 0, 8640/7877*x1, -1296/7877, 144/7877, 0; | |
| 58 -43200/43801*x1+89387/131403, 172800/43801*x1-240569/87602, -259200/43801*x1+182261/43801, 172800/43801*x1-710473/262806, -43200/43801*x1, 0, 32400/43801, -6480/43801, 720/43801 | |
| 59 ]; | |
| 60 D1(m-5:m,m-8:m) = flipud( fliplr(-D1(1:6,1:9))); | |
| 61 D1 = D1/h; | |
|
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
62 |
|
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
63 |
| 312 | 64 H = h*H; |
| 65 HI = inv(H); | |
| 66 | |
|
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
67 |
|
318
99005a80b4c2
Cleaned up d4_variable_4_min_boundary. Removed incorrect D2s from a bunch of files.
Jonatan Werpers <jonatan@werpers.com>
parents:
316
diff
changeset
|
68 D2 = []; |
|
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
69 |
| 312 | 70 % Fourth derivative, 1th order accurate at first 8 boundary points (still |
| 71 % 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
|
72 |
| 312 | 73 m4 = 7/240; |
| 74 m3 = -2/5; | |
| 75 m2 = 169/60; | |
| 76 m1 = -122/15; | |
| 77 m0 = 91/8; | |
|
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
78 |
| 312 | 79 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
|
80 |
| 312 | 81 %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
|
82 |
| 312 | 83 M4_U=[ |
| 84 0.3504379e7/0.907200e6 -0.4613983e7/0.453600e6 0.4260437e7/0.453600e6 -0.418577e6/0.113400e6 0.524579e6/0.907200e6 0.535e3/0.18144e5; | |
| 85 -0.4613983e7/0.453600e6 0.5186159e7/0.181440e6 -0.81121e5/0.2835e4 0.218845e6/0.18144e5 -0.159169e6/0.90720e5 -0.94669e5/0.907200e6; | |
| 86 0.4260437e7/0.453600e6 -0.81121e5/0.2835e4 0.147695e6/0.4536e4 -0.384457e6/0.22680e5 0.339653e6/0.90720e5 -0.18233e5/0.113400e6; | |
| 87 -0.418577e6/0.113400e6 0.218845e6/0.18144e5 -0.384457e6/0.22680e5 0.65207e5/0.4536e4 -0.22762e5/0.2835e4 0.1181753e7/0.453600e6; | |
| 88 0.524579e6/0.907200e6 -0.159169e6/0.90720e5 0.339653e6/0.90720e5 -0.22762e5/0.2835e4 0.2006171e7/0.181440e6 -0.3647647e7/0.453600e6; | |
| 89 0.535e3/0.18144e5 -0.94669e5/0.907200e6 -0.18233e5/0.113400e6 0.1181753e7/0.453600e6 -0.3647647e7/0.453600e6 0.10305271e8/0.907200e6; | |
| 90 ]; | |
| 91 | |
| 92 M4(1:6,1:6) = M4_U; | |
|
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
93 |
| 312 | 94 M4(m-5:m,m-5:m) = flipud( fliplr( M4_U ) ); |
| 95 M4 = M4/h^3; | |
|
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
96 |
| 312 | 97 D4 = HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m); |
| 98 | |
| 99 end |