Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d4_variable_6_3.m @ 322:def409c10800 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clean up of d4_variable_6_3.m
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 26 Sep 2016 09:12:08 +0200
parents 99005a80b4c2
children c0cbffcf6513
comparison
equal deleted inserted replaced
321:5c9e5ba1c1ab 322:def409c10800
18 18
19 % H?r med 7 RP ist?llet f?r 8 f?r D4 operatorn, dock samma randderivator 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 20 % Denna ?r noggrannare, och har 2a ordningens randdslutning och b?r ge 6te
21 % ordningens konvergens. Hade 2 fria parametrar att optimera 21 % ordningens konvergens. Hade 2 fria parametrar att optimera
22 22
23 % Norm
24 Hv = ones(m,1);
25 Hv(1:7) = [0.414837907e9/0.1191965760e10, 0.475278367e9/0.397321920e9, 0.13872751e8/0.12416310e8, 0.346739027e9/0.595982880e9, 0.560227469e9/0.397321920e9, 0.322971631e9/0.397321920e9, 0.616122491e9/0.595982880e9];
26 Hv(m-6:m) = rot90(Hv(1:7),2);
27 Hv = h*Hv;
28 H = spdiag(Hv, 0);
29 HI = spdiag(1./Hv, 0);
23 30
24 31
25 H=diag(ones(m,1),0); 32 % Boundary operators
26 H(1:7,1:7)=[ 33 e_l = sparse(m,1);
27 0.414837907e9/0.1191965760e10 0 0 0 0 0 0; 34 e_l(1) = 1;
28 0 0.475278367e9/0.397321920e9 0 0 0 0 0; 35 e_r = rot90(e_l, 2);
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 ];
35 36
36 H(m-6:m,m-6:m) = fliplr(flipud(H(1:7,1:7))); 37 d1_l = sparse(m,1);
38 d1_l(1:6) = [-0.137e3/0.60e2 5 -5 0.10e2/0.3e1 -0.5e1/0.4e1 0.1e1/0.5e1;]/h;
39 d1_r = -rot90(d1_l);
40
41 d2_l = sparse(m,1);
42 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;
43 d2_r = rot90(d2_l, 2);
44
45 d3_l = sparse(m,1);
46 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;
47 d3_r = -rot90(d3_l, 2);
37 48
38 49
39 e_1=zeros(m,1); 50 % Fourth derivative, 1th order accurate at first 8 boundary points
40 e_1(1)=1; 51 stencil = [7/240, -2/5, 169/60, -122/15, 91/8, -122/15, 169/60, -2/5, 7/240];
41 e_m=zeros(m,1); 52 diags = -4:4;
42 e_m(m)=1; 53 M4 = stripeMatrix(stencil, diags, m);
43
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);
49
50
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;
52 S2_1 = zeros(1,m);
53 S2_1(1:6) = S2_U;
54 S2_m = zeros(1,m);
55 S2_m(m-5:m) = fliplr(S2_U);
56
57
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;
59 S3_1 = zeros(1,m);
60 S3_1(1:6) = S3_U;
61 S3_m = zeros(1,m);
62 S3_m(m-5:m) = fliplr(-S3_U);
63
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;
68
69
70 H=h*H;
71 HI=inv(H);
72
73
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
76
77 m4 = 7/240;
78 m3 = -2/5;
79 m2 = 169/60;
80 m1 = -122/15;
81 m0 = 91/8;
82
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);
84
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));
86 54
87 M4_U = [ 55 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; 56 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; 57 -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; 58 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;
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; 61 -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; 62 -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 ]; 63 ];
96 64
97 M4(1:7,1:7) = M4_U; 65 M4(1:7,1:7) = M4_U;
66 M4(m-6:m,m-6:m) = rot90(M4_U, 2);
67 M4 = 1/h^3*M4;
98 68
99 M4(m-6:m,m-6:m) = flipud( fliplr( M4_U ) ); 69 D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r');
100 M4 = M4/h^3;
101
102 D4 = HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m);
103 end 70 end