Mercurial > repos > public > sbplib

comparison +sbp/higher_variable4_min_boundary_points.m @ 249:02423f9323c6 feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Started modifying operator files.
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 05 Sep 2016 16:38:41 +0200
parents fe26791489e0
children
comparison
equal deleted inserted replaced
248:81e0ead29431 249:02423f9323c6
1 function [H, HI, D2, D4, e_1, e_m, M4, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = higher_variable4_min_boundary_points(m,h)
1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2 %%% 4:de ordn. SBP Finita differens %%% 3 %%% 4:de ordn. SBP Finita differens %%%
3 %%% operatorer framtagna av Mark Carpenter %%% 4 %%% operatorer framtagna av Mark Carpenter %%%
4 %%% %%% 5 %%% %%%
5 %%% H (Normen) %%% 6 %%% H (Normen) %%%
9 %H?r med endast 4 randpunkter 10 %H?r med endast 4 randpunkter
10 11
11 %m=20; %problemstorlek 12 %m=20; %problemstorlek
12 %h=1/(m-1); 13 %h=1/(m-1);
13 %h=1; 14 %h=1;
14
15 % x=1:h:m*h;x=x';
16 % %c=x.^2/fac(2);
17 % x0=x.^0/fac(1);
18 % x1=x.^1/fac(1);
19 % x2=x.^2/fac(2);
20 % x3=x.^3/fac(3);
21 % x4=x.^4/fac(4);
22 % x5=x.^5/fac(5);
23 %c=ones(m,1);
24 15
25 16 Hv = ones(m,1)
26 H=diag(ones(m,1),0); 17 Hv(1:4) = [17/48 59/48 43/48 49/48];
27 H(1:4,1:4)=diag([17/48 59/48 43/48 49/48]); 18 Hv(m-3:m) = rot90(Hv(1:4),2);
28 H(m-3:m,m-3:m)=fliplr(flipud(diag([17/48 59/48 43/48 49/48]))); 19 Hv = h * Hv;
29 H=H*h; 20 H = spdiag(Hv,0);
30 HI=inv(H); 21 HI = spdiag(1./Hv,0);
31 22
32 23
24 % Boundary operators
25 e_1 = sparse(m,1);
26 e_1(1) = 1;
27 e_m = rot90(e_1, 2);
28
29 S_U=[-0.11e2 / 0.6e1 3 -0.3e1 / 0.2e1 0.1e1 / 0.3e1;]/h;
30 S_1=sparse(1,m);
31 S_1(1:4)=S_U;
32 S_m = -rot90(S_1, 2);
33
34 S2_U=[2 -5 4 -1;]/h^2;
35 S2_1=sparse(1,m);
36 S2_1(1:4)=S2_U;
37 S2_m = rot90(S2_1, 2);
38
39 S3_U=[-1 3 -3 1;]/h^3;
40 S3_1=sparse(1,m);
41 S3_1(1:4)=S3_U;
42 S3_m = rot90(S3_1, 2);
43
44
45
46
47
48
49 q_diags = [-1, 0, 1, 2];
50 q_stencil = [-1/3 -1/2 1 -1/6];
51 Qp = stripeMatrix(q_stencil, q_diags,m);
52
53 Q_U = [
54 -1/24 17/24 -1/6;
55 -13/24 -1/4 23/24;
56 1/12 -11/24 -11/24;
57 ];
58 Qp(1:3,1:3)=Q_U;
59 Qp(m-2:m,m-2:m)=rot90(Q_U,2)'; %%% This is different from standard SBP
60
61 Qm=-Qp';
62
63
33 Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2)); 64 Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2));
34 Q_U = [0 0.59e2 / 0.96e2 -0.1e1 / 0.12e2 -0.1e1 / 0.32e2; -0.59e2 / 0.96e2 0 0.59e2 / 0.96e2 0; 0.1e1 / 0.12e2 -0.59e2 / 0.96e2 0 0.59e2 / 0.96e2; 0.1e1 / 0.32e2 0 -0.59e2 / 0.96e2 0;]; 65 Q_U = [0 0.59e2 / 0.96e2 -0.1e1 / 0.12e2 -0.1e1 / 0.32e2; -0.59e2 / 0.96e2 0 0.59e2 / 0.96e2 0; 0.1e1 / 0.12e2 -0.59e2 / 0.96e2 0 0.59e2 / 0.96e2; 0.1e1 / 0.32e2 0 -0.59e2 / 0.96e2 0;];
35 Q(1:4,1:4)=Q_U; 66 Q(1:4,1:4)=Q_U;
36 Q(m-3:m,m-3:m)=flipud( fliplr(-Q_U(1:4,1:4) ) ); 67 Q(m-3:m,m-3:m)=flipud( fliplr(-Q_U(1:4,1:4) ) );
37 68
38 e_1=zeros(m,1);e_1(1)=1; 69
39 e_m=zeros(m,1);e_m(m)=1; 70
40 71
41 D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ; 72 D1=HI*(Q-1/2*e_1*e_1'+1/2*e_m*e_m') ;
42 73
43 M_U=[0.9e1 / 0.8e1 -0.59e2 / 0.48e2 0.1e1 / 0.12e2 0.1e1 / 0.48e2; -0.59e2 / 0.48e2 0.59e2 / 0.24e2 -0.59e2 / 0.48e2 0; 0.1e1 / 0.12e2 -0.59e2 / 0.48e2 0.55e2 / 0.24e2 -0.59e2 / 0.48e2; 0.1e1 / 0.48e2 0 -0.59e2 / 0.48e2 0.59e2 / 0.24e2;]; 74 M_U=[0.9e1 / 0.8e1 -0.59e2 / 0.48e2 0.1e1 / 0.12e2 0.1e1 / 0.48e2; -0.59e2 / 0.48e2 0.59e2 / 0.24e2 -0.59e2 / 0.48e2 0; 0.1e1 / 0.12e2 -0.59e2 / 0.48e2 0.55e2 / 0.24e2 -0.59e2 / 0.48e2; 0.1e1 / 0.48e2 0 -0.59e2 / 0.48e2 0.59e2 / 0.24e2;];
44 M=-(-1/12*diag(ones(m-2,1),2)+16/12*diag(ones(m-1,1),1)+16/12*diag(ones(m-1,1),-1)-1/12*diag(ones(m-2,1),-2)-30/12*diag(ones(m,1),0)); 75 M=-(-1/12*diag(ones(m-2,1),2)+16/12*diag(ones(m-1,1),1)+16/12*diag(ones(m-1,1),-1)-1/12*diag(ones(m-2,1),-2)-30/12*diag(ones(m,1),0));
45 76
46 M(1:4,1:4)=M_U; 77 M(1:4,1:4)=M_U;
47 78
48 M(m-3:m,m-3:m)=flipud( fliplr( M_U ) ); 79 M(m-3:m,m-3:m)=flipud( fliplr( M_U ) );
49 M=M/h; 80 M=M/h;
50
51 S_U=[-0.11e2 / 0.6e1 3 -0.3e1 / 0.2e1 0.1e1 / 0.3e1;]/h;
52 S_1=zeros(1,m);
53 S_1(1:4)=S_U;
54 S_m=zeros(1,m);
55 S_m(m-3:m)=fliplr(-S_U);
56
57 M=zeros(m,m); 81 M=zeros(m,m);
58 82
59 for i=4:m-3 83 for i=4:m-3
60 M(i,i-2:i+2)=[-c(i-1) / 0.6e1 + c(i-2) / 0.8e1 + c(i) / 0.8e1 -c(i-2) / 0.6e1 - c(i+1) / 0.6e1 - c(i-1) / 0.2e1 - c(i) / 0.2e1 c(i-2) / 0.24e2 + 0.5e1 / 0.6e1 * c(i-1) + 0.5e1 / 0.6e1 * c(i+1) + c(i+2) / 0.24e2 + 0.3e1 / 0.4e1 * c(i) -c(i-1) / 0.6e1 - c(i+2) / 0.6e1 - c(i) / 0.2e1 - c(i+1) / 0.2e1 -c(i+1) / 0.6e1 + c(i) / 0.8e1 + c(i+2) / 0.8e1;]; 84 M(i,i-2:i+2)=[-c(i-1) / 0.6e1 + c(i-2) / 0.8e1 + c(i) / 0.8e1 -c(i-2) / 0.6e1 - c(i+1) / 0.6e1 - c(i-1) / 0.2e1 - c(i) / 0.2e1 c(i-2) / 0.24e2 + 0.5e1 / 0.6e1 * c(i-1) + 0.5e1 / 0.6e1 * c(i+1) + c(i+2) / 0.24e2 + 0.3e1 / 0.4e1 * c(i) -c(i-1) / 0.6e1 - c(i+2) / 0.6e1 - c(i) / 0.2e1 - c(i+1) / 0.2e1 -c(i+1) / 0.6e1 + c(i) / 0.8e1 + c(i+2) / 0.8e1;];
61 end 85 end
65 M(m-5:m,m-5:m)=[c(m-7) / 0.24e2 + 0.5e1 / 0.6e1 * c(m-6) + 0.5580181e7 / 0.6692148e7 * c(m-4) + 0.4887707739997e13 / 0.119037827289528e15 * c(m-3) + 0.3e1 / 0.4e1 * c(m-5) + 0.660204843e9 / 0.13226425254392e14 * c(m-2) + 0.660204843e9 / 0.13226425254392e14 * c(m-1) -c(m-6) / 0.6e1 - 0.1618585929605e13 / 0.9919818940794e13 * c(m-3) - c(m-5) / 0.2e1 - 0.1118749e7 / 0.2230716e7 * c(m-4) - 0.13091810925e11 / 0.13226425254392e14 * c(m-2) - 0.13091810925e11 / 0.13226425254392e14 * c(m-1) -0.368395e6 / 0.2230716e7 * c(m-4) + c(m-5) / 0.8e1 + 0.48866620889e11 / 0.404890569012e12 * c(m-3) + 0.752806667e9 / 0.539854092016e12 * c(m-2) + 0.752806667e9 / 0.539854092016e12 * c(m-1) -0.3391e4 / 0.6692148e7 * c(m-4) - 0.238797444493e12 / 0.119037827289528e15 * c(m-3) + 0.33235054191e11 / 0.26452850508784e14 * c(m-2) + 0.33235054191e11 / 0.26452850508784e14 * c(m-1) -0.8673e4 / 0.2904112e7 * c(m-2) - 0.8673e4 / 0.2904112e7 * c(m-1) + 0.8673e4 / 0.1452056e7 * c(m-3) -c(m-3) / 0.392e3 + c(m-2) / 0.784e3 + c(m-1) / 0.784e3; -c(m-6) / 0.6e1 - 0.1618585929605e13 / 0.9919818940794e13 * c(m-3) - c(m-5) / 0.2e1 - 0.1118749e7 / 0.2230716e7 * c(m-4) - 0.13091810925e11 / 0.13226425254392e14 * c(m-2) - 0.13091810925e11 / 0.13226425254392e14 * c(m-1) c(m-6) / 0.24e2 + 0.5e1 / 0.6e1 * c(m-5) + 0.3896014498639e13 / 0.4959909470397e13 * c(m-3) + 0.8386761355510099813e19 / 0.128413970713633903242e21 * c(m-2) + 0.280535e6 / 0.371786e6 * c(m-4) + 0.3360696339136261875e19 / 0.171218627618178537656e21 * c(m-1) -c(m-5) / 0.6e1 - 0.4959271814984644613e19 / 0.20965546238960637264e20 * c(m-2) - 0.375177e6 / 0.743572e6 * c(m-4) - 0.13425842714e11 / 0.33740880751e11 * c(m-3) - 0.193247108773400725e18 / 0.6988515412986879088e19 * c(m-1) -0.365281640980e12 / 0.1653303156799e13 * c(m-3) + 0.564461e6 / 0.4461432e7 * c(m-4) + 0.1613976761032884305e19 / 0.7963657098519931984e19 * c(m-2) - 0.198407225513315475e18 / 0.7963657098519931984e19 * c(m-1) -0.1328188692663e13 / 0.37594290333616e14 * c(m-2) + 0.2226377963775e13 / 0.37594290333616e14 * c(m-1) - 0.8673e4 / 0.363014e6 * c(m-3) c(m-3) / 0.49e2 + 0.49579087e8 / 0.10149031312e11 * c(m-2) - 0.256702175e9 / 0.10149031312e11 * c(m-1); -0.368395e6 / 0.2230716e7 * c(m-4) + c(m-5) / 0.8e1 + 0.48866620889e11 / 0.404890569012e12 * c(m-3) + 0.752806667e9 / 0.539854092016e12 * c(m-2) + 0.752806667e9 / 0.539854092016e12 * c(m-1) -c(m-5) / 0.6e1 - 0.4959271814984644613e19 / 0.20965546238960637264e20 * c(m-2) - 0.375177e6 / 0.743572e6 * c(m-4) - 0.13425842714e11 / 0.33740880751e11 * c(m-3) - 0.193247108773400725e18 / 0.6988515412986879088e19 * c(m-1) c(m-5) / 0.24e2 + 0.1869103e7 / 0.2230716e7 * c(m-4) + 0.507284006600757858213e21 / 0.475219048083107777984e21 * c(m-2) + 0.3e1 / 0.1088e4 * c(m) + 0.31688435395e11 / 0.67481761502e11 * c(m-3) + 0.27769176016102795561e20 / 0.712828572124661666976e21 * c(m-1) -0.125059e6 / 0.743572e6 * c(m-4) + c(m) / 0.136e3 - 0.23099342648e11 / 0.101222642253e12 * c(m-3) - 0.4836340090442187227e19 / 0.5525802884687299744e19 * c(m-2) + 0.193950157930938693e18 / 0.5525802884687299744e19 * c(m-1) 0.260297319232891e15 / 0.2556411742685888e16 * c(m-2) - 0.59e2 / 0.1088e4 * c(m) - 0.106641839640553e15 / 0.1278205871342944e16 * c(m-1) + 0.26019e5 / 0.726028e6 * c(m-3) -0.1244724001e10 / 0.21126554976e11 * c(m-2) + 0.3e1 / 0.68e2 * c(m) + 0.752806667e9 / 0.21126554976e11 * c(m-1); -0.3391e4 / 0.6692148e7 * c(m-4) - 0.238797444493e12 / 0.119037827289528e15 * c(m-3) + 0.33235054191e11 / 0.26452850508784e14 * c(m-2) + 0.33235054191e11 / 0.26452850508784e14 * c(m-1) -0.365281640980e12 / 0.1653303156799e13 * c(m-3) + 0.564461e6 / 0.4461432e7 * c(m-4) + 0.1613976761032884305e19 / 0.7963657098519931984e19 * c(m-2) - 0.198407225513315475e18 / 0.7963657098519931984e19 * c(m-1) -0.125059e6 / 0.743572e6 * c(m-4) + c(m) / 0.136e3 - 0.23099342648e11 / 0.101222642253e12 * c(m-3) - 0.4836340090442187227e19 / 0.5525802884687299744e19 * c(m-2) + 0.193950157930938693e18 / 0.5525802884687299744e19 * c(m-1) 0.564461e6 / 0.13384296e8 * c(m-4) + 0.470299699916357e15 / 0.952302618316224e15 * c(m-3) + 0.550597048646198778781e21 / 0.1624586048098066124736e22 * c(m-1) + c(m) / 0.51e2 + 0.378288882302546512209e21 / 0.270764341349677687456e21 * c(m-2) -0.59e2 / 0.408e3 * c(m) - 0.29294615794607e14 / 0.29725717938208e14 * c(m-2) - 0.2234477713167e13 / 0.29725717938208e14 * c(m-1) - 0.8673e4 / 0.363014e6 * c(m-3) -0.59e2 / 0.3136e4 * c(m-3) - 0.13249937023e11 / 0.48148892736e11 * c(m-1) + 0.2e1 / 0.17e2 * c(m) + 0.2083938599e10 / 0.8024815456e10 * c(m-2); -0.8673e4 / 0.2904112e7 * c(m-2) - 0.8673e4 / 0.2904112e7 * c(m-1) + 0.8673e4 / 0.1452056e7 * c(m-3) -0.1328188692663e13 / 0.37594290333616e14 * c(m-2) + 0.2226377963775e13 / 0.37594290333616e14 * c(m-1) - 0.8673e4 / 0.363014e6 * c(m-3) 0.260297319232891e15 / 0.2556411742685888e16 * c(m-2) - 0.59e2 / 0.1088e4 * c(m) - 0.106641839640553e15 / 0.1278205871342944e16 * c(m-1) + 0.26019e5 / 0.726028e6 * c(m-3) -0.59e2 / 0.408e3 * c(m) - 0.29294615794607e14 / 0.29725717938208e14 * c(m-2) - 0.2234477713167e13 / 0.29725717938208e14 * c(m-1) - 0.8673e4 / 0.363014e6 * c(m-3) 0.9258282831623875e16 / 0.7669235228057664e16 * c(m-2) + 0.3481e4 / 0.3264e4 * c(m) + 0.228389721191751e15 / 0.1278205871342944e16 * c(m-1) + 0.8673e4 / 0.1452056e7 * c(m-3) -0.6025413881e10 / 0.21126554976e11 * c(m-2) - 0.59e2 / 0.68e2 * c(m) - 0.537416663e9 / 0.7042184992e10 * c(m-1); -c(m-3) / 0.392e3 + c(m-2) / 0.784e3 + c(m-1) / 0.784e3 c(m-3) / 0.49e2 + 0.49579087e8 / 0.10149031312e11 * c(m-2) - 0.256702175e9 / 0.10149031312e11 * c(m-1) -0.1244724001e10 / 0.21126554976e11 * c(m-2) + 0.3e1 / 0.68e2 * c(m) + 0.752806667e9 / 0.21126554976e11 * c(m-1) -0.59e2 / 0.3136e4 * c(m-3) - 0.13249937023e11 / 0.48148892736e11 * c(m-1) + 0.2e1 / 0.17e2 * c(m) + 0.2083938599e10 / 0.8024815456e10 * c(m-2) -0.6025413881e10 / 0.21126554976e11 * c(m-2) - 0.59e2 / 0.68e2 * c(m) - 0.537416663e9 / 0.7042184992e10 * c(m-1) 0.3e1 / 0.3136e4 * c(m-3) + 0.27010400129e11 / 0.345067064608e12 * c(m-2) + 0.234566387291e12 / 0.690134129216e12 * c(m-1) + 0.12e2 / 0.17e2 * c(m);]; 89 M(m-5:m,m-5:m)=[c(m-7) / 0.24e2 + 0.5e1 / 0.6e1 * c(m-6) + 0.5580181e7 / 0.6692148e7 * c(m-4) + 0.4887707739997e13 / 0.119037827289528e15 * c(m-3) + 0.3e1 / 0.4e1 * c(m-5) + 0.660204843e9 / 0.13226425254392e14 * c(m-2) + 0.660204843e9 / 0.13226425254392e14 * c(m-1) -c(m-6) / 0.6e1 - 0.1618585929605e13 / 0.9919818940794e13 * c(m-3) - c(m-5) / 0.2e1 - 0.1118749e7 / 0.2230716e7 * c(m-4) - 0.13091810925e11 / 0.13226425254392e14 * c(m-2) - 0.13091810925e11 / 0.13226425254392e14 * c(m-1) -0.368395e6 / 0.2230716e7 * c(m-4) + c(m-5) / 0.8e1 + 0.48866620889e11 / 0.404890569012e12 * c(m-3) + 0.752806667e9 / 0.539854092016e12 * c(m-2) + 0.752806667e9 / 0.539854092016e12 * c(m-1) -0.3391e4 / 0.6692148e7 * c(m-4) - 0.238797444493e12 / 0.119037827289528e15 * c(m-3) + 0.33235054191e11 / 0.26452850508784e14 * c(m-2) + 0.33235054191e11 / 0.26452850508784e14 * c(m-1) -0.8673e4 / 0.2904112e7 * c(m-2) - 0.8673e4 / 0.2904112e7 * c(m-1) + 0.8673e4 / 0.1452056e7 * c(m-3) -c(m-3) / 0.392e3 + c(m-2) / 0.784e3 + c(m-1) / 0.784e3; -c(m-6) / 0.6e1 - 0.1618585929605e13 / 0.9919818940794e13 * c(m-3) - c(m-5) / 0.2e1 - 0.1118749e7 / 0.2230716e7 * c(m-4) - 0.13091810925e11 / 0.13226425254392e14 * c(m-2) - 0.13091810925e11 / 0.13226425254392e14 * c(m-1) c(m-6) / 0.24e2 + 0.5e1 / 0.6e1 * c(m-5) + 0.3896014498639e13 / 0.4959909470397e13 * c(m-3) + 0.8386761355510099813e19 / 0.128413970713633903242e21 * c(m-2) + 0.280535e6 / 0.371786e6 * c(m-4) + 0.3360696339136261875e19 / 0.171218627618178537656e21 * c(m-1) -c(m-5) / 0.6e1 - 0.4959271814984644613e19 / 0.20965546238960637264e20 * c(m-2) - 0.375177e6 / 0.743572e6 * c(m-4) - 0.13425842714e11 / 0.33740880751e11 * c(m-3) - 0.193247108773400725e18 / 0.6988515412986879088e19 * c(m-1) -0.365281640980e12 / 0.1653303156799e13 * c(m-3) + 0.564461e6 / 0.4461432e7 * c(m-4) + 0.1613976761032884305e19 / 0.7963657098519931984e19 * c(m-2) - 0.198407225513315475e18 / 0.7963657098519931984e19 * c(m-1) -0.1328188692663e13 / 0.37594290333616e14 * c(m-2) + 0.2226377963775e13 / 0.37594290333616e14 * c(m-1) - 0.8673e4 / 0.363014e6 * c(m-3) c(m-3) / 0.49e2 + 0.49579087e8 / 0.10149031312e11 * c(m-2) - 0.256702175e9 / 0.10149031312e11 * c(m-1); -0.368395e6 / 0.2230716e7 * c(m-4) + c(m-5) / 0.8e1 + 0.48866620889e11 / 0.404890569012e12 * c(m-3) + 0.752806667e9 / 0.539854092016e12 * c(m-2) + 0.752806667e9 / 0.539854092016e12 * c(m-1) -c(m-5) / 0.6e1 - 0.4959271814984644613e19 / 0.20965546238960637264e20 * c(m-2) - 0.375177e6 / 0.743572e6 * c(m-4) - 0.13425842714e11 / 0.33740880751e11 * c(m-3) - 0.193247108773400725e18 / 0.6988515412986879088e19 * c(m-1) c(m-5) / 0.24e2 + 0.1869103e7 / 0.2230716e7 * c(m-4) + 0.507284006600757858213e21 / 0.475219048083107777984e21 * c(m-2) + 0.3e1 / 0.1088e4 * c(m) + 0.31688435395e11 / 0.67481761502e11 * c(m-3) + 0.27769176016102795561e20 / 0.712828572124661666976e21 * c(m-1) -0.125059e6 / 0.743572e6 * c(m-4) + c(m) / 0.136e3 - 0.23099342648e11 / 0.101222642253e12 * c(m-3) - 0.4836340090442187227e19 / 0.5525802884687299744e19 * c(m-2) + 0.193950157930938693e18 / 0.5525802884687299744e19 * c(m-1) 0.260297319232891e15 / 0.2556411742685888e16 * c(m-2) - 0.59e2 / 0.1088e4 * c(m) - 0.106641839640553e15 / 0.1278205871342944e16 * c(m-1) + 0.26019e5 / 0.726028e6 * c(m-3) -0.1244724001e10 / 0.21126554976e11 * c(m-2) + 0.3e1 / 0.68e2 * c(m) + 0.752806667e9 / 0.21126554976e11 * c(m-1); -0.3391e4 / 0.6692148e7 * c(m-4) - 0.238797444493e12 / 0.119037827289528e15 * c(m-3) + 0.33235054191e11 / 0.26452850508784e14 * c(m-2) + 0.33235054191e11 / 0.26452850508784e14 * c(m-1) -0.365281640980e12 / 0.1653303156799e13 * c(m-3) + 0.564461e6 / 0.4461432e7 * c(m-4) + 0.1613976761032884305e19 / 0.7963657098519931984e19 * c(m-2) - 0.198407225513315475e18 / 0.7963657098519931984e19 * c(m-1) -0.125059e6 / 0.743572e6 * c(m-4) + c(m) / 0.136e3 - 0.23099342648e11 / 0.101222642253e12 * c(m-3) - 0.4836340090442187227e19 / 0.5525802884687299744e19 * c(m-2) + 0.193950157930938693e18 / 0.5525802884687299744e19 * c(m-1) 0.564461e6 / 0.13384296e8 * c(m-4) + 0.470299699916357e15 / 0.952302618316224e15 * c(m-3) + 0.550597048646198778781e21 / 0.1624586048098066124736e22 * c(m-1) + c(m) / 0.51e2 + 0.378288882302546512209e21 / 0.270764341349677687456e21 * c(m-2) -0.59e2 / 0.408e3 * c(m) - 0.29294615794607e14 / 0.29725717938208e14 * c(m-2) - 0.2234477713167e13 / 0.29725717938208e14 * c(m-1) - 0.8673e4 / 0.363014e6 * c(m-3) -0.59e2 / 0.3136e4 * c(m-3) - 0.13249937023e11 / 0.48148892736e11 * c(m-1) + 0.2e1 / 0.17e2 * c(m) + 0.2083938599e10 / 0.8024815456e10 * c(m-2); -0.8673e4 / 0.2904112e7 * c(m-2) - 0.8673e4 / 0.2904112e7 * c(m-1) + 0.8673e4 / 0.1452056e7 * c(m-3) -0.1328188692663e13 / 0.37594290333616e14 * c(m-2) + 0.2226377963775e13 / 0.37594290333616e14 * c(m-1) - 0.8673e4 / 0.363014e6 * c(m-3) 0.260297319232891e15 / 0.2556411742685888e16 * c(m-2) - 0.59e2 / 0.1088e4 * c(m) - 0.106641839640553e15 / 0.1278205871342944e16 * c(m-1) + 0.26019e5 / 0.726028e6 * c(m-3) -0.59e2 / 0.408e3 * c(m) - 0.29294615794607e14 / 0.29725717938208e14 * c(m-2) - 0.2234477713167e13 / 0.29725717938208e14 * c(m-1) - 0.8673e4 / 0.363014e6 * c(m-3) 0.9258282831623875e16 / 0.7669235228057664e16 * c(m-2) + 0.3481e4 / 0.3264e4 * c(m) + 0.228389721191751e15 / 0.1278205871342944e16 * c(m-1) + 0.8673e4 / 0.1452056e7 * c(m-3) -0.6025413881e10 / 0.21126554976e11 * c(m-2) - 0.59e2 / 0.68e2 * c(m) - 0.537416663e9 / 0.7042184992e10 * c(m-1); -c(m-3) / 0.392e3 + c(m-2) / 0.784e3 + c(m-1) / 0.784e3 c(m-3) / 0.49e2 + 0.49579087e8 / 0.10149031312e11 * c(m-2) - 0.256702175e9 / 0.10149031312e11 * c(m-1) -0.1244724001e10 / 0.21126554976e11 * c(m-2) + 0.3e1 / 0.68e2 * c(m) + 0.752806667e9 / 0.21126554976e11 * c(m-1) -0.59e2 / 0.3136e4 * c(m-3) - 0.13249937023e11 / 0.48148892736e11 * c(m-1) + 0.2e1 / 0.17e2 * c(m) + 0.2083938599e10 / 0.8024815456e10 * c(m-2) -0.6025413881e10 / 0.21126554976e11 * c(m-2) - 0.59e2 / 0.68e2 * c(m) - 0.537416663e9 / 0.7042184992e10 * c(m-1) 0.3e1 / 0.3136e4 * c(m-3) + 0.27010400129e11 / 0.345067064608e12 * c(m-2) + 0.234566387291e12 / 0.690134129216e12 * c(m-1) + 0.12e2 / 0.17e2 * c(m);];
66 90
67 M=M/h; 91 M=M/h;
68 92
69 D2=HI*(-M-c(1)*e_1*S_1+c(m)*e_m*S_m); 93 D2=HI*(-M-c(1)*e_1*S_1+c(m)*e_m*S_m);
70
71 S2_U=[2 -5 4 -1;]/h^2;
72 S2_1=zeros(1,m);
73 S2_1(1:4)=S2_U;
74 S2_m=zeros(1,m);
75 S2_m(m-3:m)=fliplr(S2_U);
76
77 m3=-1/6;m2=2;m1=-13/2;m0=28/3; 94 m3=-1/6;m2=2;m1=-13/2;m0=28/3;
78 M4=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); 95 M4=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);
79 96
80 %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));
81 97
82 M4_U=[0.8e1 / 0.3e1 -0.37e2 / 0.6e1 0.13e2 / 0.3e1 -0.5e1 / 0.6e1; -0.37e2 / 0.6e1 0.47e2 / 0.3e1 -13 0.11e2 / 0.3e1; 0.13e2 / 0.3e1 -13 0.44e2 / 0.3e1 -0.47e2 / 0.6e1; -0.5e1 / 0.6e1 0.11e2 / 0.3e1 -0.47e2 / 0.6e1 0.29e2 / 0.3e1;]; 98 M4_U=[0.8e1 / 0.3e1 -0.37e2 / 0.6e1 0.13e2 / 0.3e1 -0.5e1 / 0.6e1; -0.37e2 / 0.6e1 0.47e2 / 0.3e1 -13 0.11e2 / 0.3e1; 0.13e2 / 0.3e1 -13 0.44e2 / 0.3e1 -0.47e2 / 0.6e1; -0.5e1 / 0.6e1 0.11e2 / 0.3e1 -0.47e2 / 0.6e1 0.29e2 / 0.3e1;];
83 M4(1:4,1:4)=M4_U; 99 M4(1:4,1:4)=M4_U;
84 100
85 M4(m-3:m,m-3:m)=flipud( fliplr( M4_U ) ); 101 M4(m-3:m,m-3:m)=flipud( fliplr( M4_U ) );
86 M4=M4/h^3; 102 M4=M4/h^3;
87
88 S3_U=[-1 3 -3 1;]/h^3;
89 S3_1=zeros(1,m);
90 S3_1(1:4)=S3_U;
91 S3_m=zeros(1,m);
92 S3_m(m-3:m)=fliplr(-S3_U);
93 103
94 D4=HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m); 104 D4=HI*(M4-e_1*S3_1+e_m*S3_m + S_1'*S2_1-S_m'*S2_m);
95
96 % % Test
97 % d3=[-1 3 -3 1];
98 % DD_3=(-diag(ones(m-2,1),-2)+3*diag(ones(m-1,1),-1)-3*diag(ones(m,1),0)+ ...
99 % diag(ones(m-1,1),1));
100 % DD_3(1:2,1:4)=[d3;d3];
101 % DD_3(m,m-3:m)=[d3];
102 %
103 % nnn=round(0.1/h);
104 % x0=h;
105 % x1=nnn*h;
106 % d=h; % Start level
107 % g=1; % interrior level
108 %
109 % c0 = (-10*x1^3*x0^2*d+5*x1^4*x0*d-d*x1^5+10*x1^2*x0^3*g-5*x1*x0^4*g+x0^5* ...
110 % g)/(10*x1^2*x0^3-10*x0^2*x1^3+5*x0*x1^4-5*x1*x0^4+x0^5-x1^5);
111 % c1 = -30*(-d+g)*x1^2*x0^2/(10*x1^2*x0^3-10*x0^2*x1^3+5*x0*x1^4-5*x1*x0^4+x0^5-x1^5);
112 % c2 = 30*(-d+g)*x1*x0*(x0+x1)/(10*x1^2*x0^3-10*x0^2*x1^3+5*x0*x1^4-5*x1* ...
113 % x0^4+x0^5-x1^5);
114 % c3 = -10*(-d+g)*(x0^2+4*x1*x0+x1^2)/(10*x1^2*x0^3-10*x0^2*x1^3+5*x0*x1^4- ...
115 % 5*x1*x0^4+x0^5-x1^5);
116 % c4 = 15*(-d+g)*(x0+x1)/(10*x1^2*x0^3-10*x0^2*x1^3+5*x0*x1^4-5*x1*x0^4+ ...
117 % x0^5-x1^5);
118 % c5 = -6*(-d+g)/(10*x1^2*x0^3-10*x0^2*x1^3+5*x0*x1^4-5*x1*x0^4+x0^5-x1^5);
119 %
120 % xt=x0:h:x1;
121 % ct=polyval([c5 c4 c3 c2 c1 c0],xt);
122 %
123 % BB1=zeros(m);
124 % BB1(1:nnn,1:nnn)=diag(ct);
125 % BB1(m-nnn+1:m,m-nnn+1:m)=diag(fliplr(ct));
126 % BB1(nnn+1:m-nnn,nnn+1:m-nnn)=g*eye(length(nnn+1:m-nnn));
127 %
128 % DI=(1/h^3)*DD_3'*BB1*DD_3;
129 %
130 % I=eye(m);
131 % MM1=kron(H,D2'*H*D2)+kron(D2'*H*D2,H)+kron(D2'*H,H*D2)+kron(H*D2,D2'*H);
132 % ee1=eig(MM1);min(real(ee1)),max(real(ee1))
133 %
134 %
135 % MM2=kron(H,M4)+kron(M4,H)+kron(D2'*H,H*D2)+kron(H*D2,D2'*H);
136 % MM3=kron(H,M4+DI)+kron(M4+DI,H)+kron(D2'*H,H*D2)+kron(H*D2,D2'*H);
137 % ee3=eig(MM3);min(real(ee3)),max(real(ee3))
138