Mercurial > repos > public > sbplib
comparison +sbp/+implementations/d4_lonely_6_min_boundary_points.m @ 325:72468bc9b63f feature/beams
Renamed some operator implementations.
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Mon, 26 Sep 2016 09:55:16 +0200 |
| parents | +sbp/+implementations/d4_variable_6_min_boundary_points.m@c0cbffcf6513 |
| children | b19e142fcae1 |
comparison
equal
deleted
inserted
replaced
| 324:c0cbffcf6513 | 325:72468bc9b63f |
|---|---|
| 1 function [H, HI, 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% | |
| 18 | |
| 19 % H?r med 6 RP ist?llet f?r 8 f?r D4 operatorn, dock samma randderivator | |
| 20 | |
| 21 BP = 6; | |
| 22 if(m<2*BP) | |
| 23 error(['Operator requires at least ' num2str(2*BP) ' grid points']); | |
| 24 end | |
| 25 | |
| 26 % Norm | |
| 27 Hv = ones(m,1); | |
| 28 Hv(1:6) = [13649/43200,12013/8640,2711/4320,5359/4320,7877/8640, 43801/43200]; | |
| 29 Hv(m-5:m) = rot90(Hv(1:6),2); | |
| 30 Hv = h*Hv; | |
| 31 H = spdiag(Hv, 0); | |
| 32 HI = spdiag(1./Hv, 0); | |
| 33 | |
| 34 | |
| 35 % Boundary operators | |
| 36 e_l = sparse(m,1); | |
| 37 e_l(1) = 1; | |
| 38 e_r = rot90(e_l, 2); | |
| 39 | |
| 40 d1_l = sparse(m,1); | |
| 41 d1_l(1:5) = [-25/12, 4, -3, 4/3, -1/4]/h; | |
| 42 d1_r = -rot90(d1_l); | |
| 43 | |
| 44 d2_l = sparse(m,1); | |
| 45 d2_l(1:5) = [0.35e2/0.12e2 -0.26e2/0.3e1 0.19e2/0.2e1 -0.14e2/0.3e1 0.11e2/0.12e2;]/h^2; | |
| 46 d2_r = rot90(d2_l, 2); | |
| 47 | |
| 48 d3_l = sparse(m,1); | |
| 49 d3_l(1:5) = [-0.5e1/0.2e1 9 -12 7 -0.3e1/0.2e1;]/h^3; | |
| 50 d3_r = -rot90(d3_l, 2); | |
| 51 | |
| 52 | |
| 53 % Fourth derivative, 1th order accurate at first 8 boundary points (still | |
| 54 % yield 5th order convergence if stable: for example u_tt=-u_xxxx | |
| 55 | |
| 56 stencil = [7/240, -2/5, 169/60, -122/15, 91/8, -122/15, 169/60, -2/5, 7/240]; | |
| 57 diags = -4:4; | |
| 58 M4 = stripeMatrix(stencil, diags, m); | |
| 59 | |
| 60 M4_U=[ | |
| 61 0.3504379e7/0.907200e6 -0.4613983e7/0.453600e6 0.4260437e7/0.453600e6 -0.418577e6/0.113400e6 0.524579e6/0.907200e6 0.535e3/0.18144e5; | |
| 62 -0.4613983e7/0.453600e6 0.5186159e7/0.181440e6 -0.81121e5/0.2835e4 0.218845e6/0.18144e5 -0.159169e6/0.90720e5 -0.94669e5/0.907200e6; | |
| 63 0.4260437e7/0.453600e6 -0.81121e5/0.2835e4 0.147695e6/0.4536e4 -0.384457e6/0.22680e5 0.339653e6/0.90720e5 -0.18233e5/0.113400e6; | |
| 64 -0.418577e6/0.113400e6 0.218845e6/0.18144e5 -0.384457e6/0.22680e5 0.65207e5/0.4536e4 -0.22762e5/0.2835e4 0.1181753e7/0.453600e6; | |
| 65 0.524579e6/0.907200e6 -0.159169e6/0.90720e5 0.339653e6/0.90720e5 -0.22762e5/0.2835e4 0.2006171e7/0.181440e6 -0.3647647e7/0.453600e6; | |
| 66 0.535e3/0.18144e5 -0.94669e5/0.907200e6 -0.18233e5/0.113400e6 0.1181753e7/0.453600e6 -0.3647647e7/0.453600e6 0.10305271e8/0.907200e6; | |
| 67 ]; | |
| 68 | |
| 69 M4(1:6,1:6) = M4_U; | |
| 70 M4(m-5:m,m-5:m) = rot90(M4_U, 2); | |
| 71 M4 = 1/h^3*M4; | |
| 72 | |
| 73 D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r'); | |
| 74 end |