Mercurial > repos > public > sbplib
annotate +sbp/+implementations/d4_variable_6_2.m @ 321:5c9e5ba1c1ab feature/beams
Clean up of d4_variable_6_2.m
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Mon, 26 Sep 2016 09:04:57 +0200 |
parents | 99005a80b4c2 |
children | c0cbffcf6513 |
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_2(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 |
20 % Denna ?r noggrannare, och har 2a ordningens randdslutning och b?r ge 6te | |
21 % ordningens konvergens. Hade dock ingen fri parameter att optimera | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
22 |
321
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
23 BP = 6; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
24 if(m<2*BP) |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
25 error(['Operator requires at least ' num2str(2*BP) ' grid points']); |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
26 end |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
27 |
321
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
28 % Norm |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
29 Hv = ones(m,1); |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
30 Hv(1:6) = [0.181e3/0.576e3, 0.1343e4/0.960e3, 0.293e3/0.480e3, 0.1811e4/0.1440e4, 0.289e3/0.320e3, 0.65e2/0.64e2]; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
31 Hv(m-5:m) = rot90(Hv(1:6),2); |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
32 Hv = h*Hv; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
33 H = spdiag(Hv, 0); |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
34 HI = spdiag(1./Hv, 0); |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
35 |
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
36 |
321
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
37 % Boundary operators |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
38 e_l = sparse(m,1); |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
39 e_l(1) = 1; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
40 e_r = rot90(e_l, 2); |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
41 |
321
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
42 d1_l = sparse(m,1); |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
43 d1_l(1:6) = [-0.137e3/0.60e2 5 -5 0.10e2/0.3e1 -0.5e1/0.4e1 0.1e1/0.5e1;]/h; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
44 d1_r = -rot90(d1_l); |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
45 |
321
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
46 d2_l = sparse(m,1); |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
47 d2_l(1:6) = [0.15e2/0.4e1 -0.77e2/0.6e1 0.107e3/0.6e1 -13 0.61e2/0.12e2 -0.5e1/0.6e1;]/h^2; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
48 d2_r = rot90(d2_l, 2); |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
49 |
321
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
50 d3_l = sparse(m,1); |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
51 d3_l(1:6) = [-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; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
52 d3_r = -rot90(d3_l, 2); |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
53 |
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
54 |
312 | 55 % Fourth derivative, 1th order accurate at first 8 boundary points (still |
56 % yield 5th order convergence if stable: for example u_tt = -u_xxxx | |
321
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
57 stencil = [7/240, -2/5, 169/60, -122/15, 91/8, -122/15, 169/60, -2/5, 7/240]; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
58 diags = -4:4; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
59 M4 = stripeMatrix(stencil, diags, m); |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
60 |
312 | 61 M4_U = [ |
62 0.1009e4/0.192e3 -0.7657e4/0.480e3 0.9307e4/0.480e3 -0.509e3/0.40e2 0.4621e4/0.960e3 -0.25e2/0.32e2; | |
63 -0.7657e4/0.480e3 0.49513e5/0.960e3 -0.4007e4/0.60e2 0.21799e5/0.480e3 -0.8171e4/0.480e3 0.2657e4/0.960e3; | |
64 0.9307e4/0.480e3 -0.4007e4/0.60e2 0.1399e4/0.15e2 -0.2721e4/0.40e2 0.12703e5/0.480e3 -0.521e3/0.120e3; | |
65 -0.509e3/0.40e2 0.21799e5/0.480e3 -0.2721e4/0.40e2 0.3349e4/0.60e2 -0.389e3/0.15e2 0.559e3/0.96e2; | |
66 0.4621e4/0.960e3 -0.8171e4/0.480e3 0.12703e5/0.480e3 -0.389e3/0.15e2 0.17857e5/0.960e3 -0.1499e4/0.160e3; | |
67 -0.25e2/0.32e2 0.2657e4/0.960e3 -0.521e3/0.120e3 0.559e3/0.96e2 -0.1499e4/0.160e3 0.2225e4/0.192e3; | |
68 ]; | |
69 | |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
70 |
321
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
71 M4(1:6,1:6) = M4_U; |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
72 M4(m-5:m,m-5:m) = rot90(M4_U, 2); |
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
73 M4 = 1/h^3*M4; |
246
fe26791489e0
Added a bunch of new operators. Still non-functional.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff
changeset
|
74 |
321
5c9e5ba1c1ab
Clean up of d4_variable_6_2.m
Jonatan Werpers <jonatan@werpers.com>
parents:
318
diff
changeset
|
75 D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r'); |
312 | 76 end |