Mercurial > repos > public > sbplib
diff +sbp/+implementations/d1_noneq_12.m @ 261:6009f2712d13 operator_remake
Moved and renamned all implementations.
author | Martin Almquist <martin.almquist@it.uu.se> |
---|---|
date | Thu, 08 Sep 2016 15:35:45 +0200 |
parents | |
children | bfa130b7abf6 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d1_noneq_12.m Thu Sep 08 15:35:45 2016 +0200 @@ -0,0 +1,354 @@ +function [D1,H,x,h] = d1_noneq_12(N,L) + +% L: Domain length +% N: Number of grid points +if(nargin < 2) + L = 1; +end + +% BP: Number of boundary points +% m: Number of nonequidistant spacings +% order: Accuracy of interior stencil +BP = 12; +m = 6; +order = 12; + +%%%% Non-equidistant grid points %%%%% +x0 = 0.0000000000000e+00; +x1 = 3.6098032343909e-01; +x2 = 1.1634317168086e+00; +x3 = 2.2975905356987e+00; +x4 = 3.6057529790929e+00; +x5 = 4.8918275675510e+00; +x6 = 6.0000000000000e+00; +x7 = 7.0000000000000e+00; +x8 = 8.0000000000000e+00; +x9 = 9.0000000000000e+00; +x10 = 1.0000000000000e+01; +x11 = 1.1000000000000e+01; +x12 = 1.2000000000000e+01; + +xb = zeros(m+1,1); +for i = 0:m + xb(i+1) = eval(['x' num2str(i)]); +end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%% Compute h %%%%%%%%%% +h = L/(2*xb(end) + N-1-2*m); +%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%% Define grid %%%%%%%% +x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; +%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%% Norm matrix %%%%%%%% +P = zeros(BP,1); +%#ok<*NASGU> +P0 = 1.0000000000011e-01; +P1 = 5.9616216757547e-01; +P2 = 9.9065699844442e-01; +P3 = 1.2512548713913e+00; +P4 = 1.3316678989403e+00; +P5 = 1.2093375037721e+00; +P6 = 1.0236491595704e+00; +P7 = 9.9685258909811e-01; +P8 = 1.0004766563923e+00; +P9 = 9.9993617879146e-01; +P10 = 1.0000063122914e+00; +P11 = 9.9999966373260e-01; + +for i = 0:BP-1 + P(i+1) = eval(['P' num2str(i)]); +end + +H = ones(N,1); +H(1:BP) = P; +H(end-BP+1:end) = flip(P); +H = spdiags(h*H,0,N,N); +%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%% Q matrix %%%%%%%%%%% + +% interior stencil +switch order + case 2 + d = [-1/2,0,1/2]; + case 4 + d = [1/12,-2/3,0,2/3,-1/12]; + case 6 + d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; + case 8 + d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; + case 10 + d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; + case 12 + d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; +end +d = repmat(d,N,1); +Q = spdiags(d,-order/2:order/2,N,N); + +% Boundaries +Q0_0 = -5.0000000000000e-01; +Q0_1 = 6.7597560728423e-01; +Q0_2 = -2.6859785384416e-01; +Q0_3 = 1.4850302678903e-01; +Q0_4 = -8.7976689586154e-02; +Q0_5 = 4.1833336322613e-02; +Q0_6 = -2.2216684976993e-03; +Q0_7 = -1.5910034062022e-02; +Q0_8 = 1.1296706376589e-02; +Q0_9 = -3.1823678285130e-03; +Q0_10 = 2.4843594063649e-04; +Q0_11 = 3.1501105449828e-05; +Q0_12 = 0.0000000000000e+00; +Q0_13 = 0.0000000000000e+00; +Q0_14 = 0.0000000000000e+00; +Q0_15 = 0.0000000000000e+00; +Q0_16 = 0.0000000000000e+00; +Q0_17 = 0.0000000000000e+00; +Q1_0 = -6.7597560728423e-01; +Q1_1 = 0.0000000000000e+00; +Q1_2 = 9.5424013647146e-01; +Q1_3 = -4.3389334603464e-01; +Q1_4 = 2.4285669347653e-01; +Q1_5 = -1.1443465137214e-01; +Q1_6 = 8.5942765682435e-03; +Q1_7 = 4.0290424215772e-02; +Q1_8 = -2.9396383714543e-02; +Q1_9 = 8.5601827834256e-03; +Q1_10 = -7.8128092862319e-04; +Q1_11 = -6.0444181254875e-05; +Q1_12 = 0.0000000000000e+00; +Q1_13 = 0.0000000000000e+00; +Q1_14 = 0.0000000000000e+00; +Q1_15 = 0.0000000000000e+00; +Q1_16 = 0.0000000000000e+00; +Q1_17 = 0.0000000000000e+00; +Q2_0 = 2.6859785384416e-01; +Q2_1 = -9.5424013647146e-01; +Q2_2 = 0.0000000000000e+00; +Q2_3 = 9.7065114311923e-01; +Q2_4 = -4.3205328628292e-01; +Q2_5 = 1.9549970932735e-01; +Q2_6 = -2.4406885385172e-02; +Q2_7 = -5.5737279079895e-02; +Q2_8 = 4.3772338637753e-02; +Q2_9 = -1.3727655130726e-02; +Q2_10 = 1.6271304373071e-03; +Q2_11 = 1.7066984372933e-05; +Q2_12 = 0.0000000000000e+00; +Q2_13 = 0.0000000000000e+00; +Q2_14 = 0.0000000000000e+00; +Q2_15 = 0.0000000000000e+00; +Q2_16 = 0.0000000000000e+00; +Q2_17 = 0.0000000000000e+00; +Q3_0 = -1.4850302678903e-01; +Q3_1 = 4.3389334603464e-01; +Q3_2 = -9.7065114311923e-01; +Q3_3 = 0.0000000000000e+00; +Q3_4 = 9.5375878629204e-01; +Q3_5 = -3.6113954384951e-01; +Q3_6 = 6.9749289223875e-02; +Q3_7 = 6.5161366516465e-02; +Q3_8 = -6.0325702283960e-02; +Q3_9 = 2.1188913621662e-02; +Q3_10 = -3.2632650250470e-03; +Q3_11 = 1.3097937809499e-04; +Q3_12 = 0.0000000000000e+00; +Q3_13 = 0.0000000000000e+00; +Q3_14 = 0.0000000000000e+00; +Q3_15 = 0.0000000000000e+00; +Q3_16 = 0.0000000000000e+00; +Q3_17 = 0.0000000000000e+00; +Q4_0 = 8.7976689586154e-02; +Q4_1 = -2.4285669347653e-01; +Q4_2 = 4.3205328628292e-01; +Q4_3 = -9.5375878629204e-01; +Q4_4 = 0.0000000000000e+00; +Q4_5 = 8.8676146394834e-01; +Q4_6 = -2.1292503103800e-01; +Q4_7 = -4.6037018833218e-02; +Q4_8 = 7.4338719466734e-02; +Q4_9 = -3.1217656663809e-02; +Q4_10 = 6.1239492854797e-03; +Q4_11 = -4.5892226603067e-04; +Q4_12 = 0.0000000000000e+00; +Q4_13 = 0.0000000000000e+00; +Q4_14 = 0.0000000000000e+00; +Q4_15 = 0.0000000000000e+00; +Q4_16 = 0.0000000000000e+00; +Q4_17 = 0.0000000000000e+00; +Q5_0 = -4.1833336322613e-02; +Q5_1 = 1.1443465137214e-01; +Q5_2 = -1.9549970932735e-01; +Q5_3 = 3.6113954384951e-01; +Q5_4 = -8.8676146394834e-01; +Q5_5 = 0.0000000000000e+00; +Q5_6 = 7.7461223007026e-01; +Q5_7 = -1.0609547334165e-01; +Q5_8 = -4.4853791547749e-02; +Q5_9 = 3.2436468405486e-02; +Q5_10 = -8.4387621360184e-03; +Q5_11 = 8.5964292632428e-04; +Q5_12 = 0.0000000000000e+00; +Q5_13 = 0.0000000000000e+00; +Q5_14 = 0.0000000000000e+00; +Q5_15 = 0.0000000000000e+00; +Q5_16 = 0.0000000000000e+00; +Q5_17 = 0.0000000000000e+00; +Q6_0 = 2.2216684976993e-03; +Q6_1 = -8.5942765682435e-03; +Q6_2 = 2.4406885385172e-02; +Q6_3 = -6.9749289223875e-02; +Q6_4 = 2.1292503103800e-01; +Q6_5 = -7.7461223007026e-01; +Q6_6 = 0.0000000000000e+00; +Q6_7 = 7.4758103262966e-01; +Q6_8 = -1.5730779067906e-01; +Q6_9 = 2.6517620342970e-02; +Q6_10 = -4.3175367549700e-03; +Q6_11 = 1.1092605832824e-03; +Q6_12 = -1.8037518037522e-04; +Q6_13 = 0.0000000000000e+00; +Q6_14 = 0.0000000000000e+00; +Q6_15 = 0.0000000000000e+00; +Q6_16 = 0.0000000000000e+00; +Q6_17 = 0.0000000000000e+00; +Q7_0 = 1.5910034062022e-02; +Q7_1 = -4.0290424215772e-02; +Q7_2 = 5.5737279079895e-02; +Q7_3 = -6.5161366516465e-02; +Q7_4 = 4.6037018833218e-02; +Q7_5 = 1.0609547334165e-01; +Q7_6 = -7.4758103262966e-01; +Q7_7 = 0.0000000000000e+00; +Q7_8 = 8.0975719267918e-01; +Q7_9 = -2.3568822398349e-01; +Q7_10 = 6.9373143801571e-02; +Q7_11 = -1.6606121869177e-02; +Q7_12 = 2.5974025974031e-03; +Q7_13 = -1.8037518037522e-04; +Q7_14 = 0.0000000000000e+00; +Q7_15 = 0.0000000000000e+00; +Q7_16 = 0.0000000000000e+00; +Q7_17 = 0.0000000000000e+00; +Q8_0 = -1.1296706376589e-02; +Q8_1 = 2.9396383714543e-02; +Q8_2 = -4.3772338637753e-02; +Q8_3 = 6.0325702283960e-02; +Q8_4 = -7.4338719466734e-02; +Q8_5 = 4.4853791547749e-02; +Q8_6 = 1.5730779067906e-01; +Q8_7 = -8.0975719267918e-01; +Q8_8 = 0.0000000000000e+00; +Q8_9 = 8.4765775072084e-01; +Q8_10 = -2.6369594097148e-01; +Q8_11 = 7.8759594625702e-02; +Q8_12 = -1.7857142857146e-02; +Q8_13 = 2.5974025974031e-03; +Q8_14 = -1.8037518037522e-04; +Q8_15 = 0.0000000000000e+00; +Q8_16 = 0.0000000000000e+00; +Q8_17 = 0.0000000000000e+00; +Q9_0 = 3.1823678285130e-03; +Q9_1 = -8.5601827834256e-03; +Q9_2 = 1.3727655130726e-02; +Q9_3 = -2.1188913621662e-02; +Q9_4 = 3.1217656663809e-02; +Q9_5 = -3.2436468405486e-02; +Q9_6 = -2.6517620342970e-02; +Q9_7 = 2.3568822398349e-01; +Q9_8 = -8.4765775072084e-01; +Q9_9 = 0.0000000000000e+00; +Q9_10 = 8.5631774953989e-01; +Q9_11 = -2.6769768119702e-01; +Q9_12 = 7.9365079365093e-02; +Q9_13 = -1.7857142857146e-02; +Q9_14 = 2.5974025974031e-03; +Q9_15 = -1.8037518037522e-04; +Q9_16 = 0.0000000000000e+00; +Q9_17 = 0.0000000000000e+00; +Q10_0 = -2.4843594063649e-04; +Q10_0 = -2.4843594063649e-04; +Q10_1 = 7.8128092862319e-04; +Q10_1 = 7.8128092862319e-04; +Q10_2 = -1.6271304373071e-03; +Q10_2 = -1.6271304373071e-03; +Q10_3 = 3.2632650250470e-03; +Q10_3 = 3.2632650250470e-03; +Q10_4 = -6.1239492854797e-03; +Q10_4 = -6.1239492854797e-03; +Q10_5 = 8.4387621360184e-03; +Q10_5 = 8.4387621360184e-03; +Q10_6 = 4.3175367549700e-03; +Q10_6 = 4.3175367549700e-03; +Q10_7 = -6.9373143801571e-02; +Q10_7 = -6.9373143801571e-02; +Q10_8 = 2.6369594097148e-01; +Q10_8 = 2.6369594097148e-01; +Q10_9 = -8.5631774953989e-01; +Q10_9 = -8.5631774953989e-01; +Q10_10 = 0.0000000000000e+00; +Q10_10 = 0.0000000000000e+00; +Q10_11 = 8.5712580212095e-01; +Q10_11 = 8.5712580212095e-01; +Q10_12 = -2.6785714285718e-01; +Q10_12 = -2.6785714285718e-01; +Q10_13 = 7.9365079365093e-02; +Q10_13 = 7.9365079365093e-02; +Q10_14 = -1.7857142857146e-02; +Q10_14 = -1.7857142857146e-02; +Q10_15 = 2.5974025974031e-03; +Q10_15 = 2.5974025974031e-03; +Q10_16 = -1.8037518037522e-04; +Q10_16 = -1.8037518037522e-04; +Q10_17 = 0.0000000000000e+00; +Q10_17 = 0.0000000000000e+00; +Q11_0 = -3.1501105449828e-05; +Q11_0 = -3.1501105449828e-05; +Q11_1 = 6.0444181254875e-05; +Q11_1 = 6.0444181254875e-05; +Q11_2 = -1.7066984372933e-05; +Q11_2 = -1.7066984372933e-05; +Q11_3 = -1.3097937809499e-04; +Q11_3 = -1.3097937809499e-04; +Q11_4 = 4.5892226603067e-04; +Q11_4 = 4.5892226603067e-04; +Q11_5 = -8.5964292632428e-04; +Q11_5 = -8.5964292632428e-04; +Q11_6 = -1.1092605832824e-03; +Q11_6 = -1.1092605832824e-03; +Q11_7 = 1.6606121869177e-02; +Q11_7 = 1.6606121869177e-02; +Q11_8 = -7.8759594625702e-02; +Q11_8 = -7.8759594625702e-02; +Q11_9 = 2.6769768119702e-01; +Q11_9 = 2.6769768119702e-01; +Q11_10 = -8.5712580212095e-01; +Q11_10 = -8.5712580212095e-01; +Q11_11 = 0.0000000000000e+00; +Q11_11 = 0.0000000000000e+00; +Q11_12 = 8.5714285714289e-01; +Q11_12 = 8.5714285714289e-01; +Q11_13 = -2.6785714285718e-01; +Q11_13 = -2.6785714285718e-01; +Q11_14 = 7.9365079365093e-02; +Q11_14 = 7.9365079365093e-02; +Q11_15 = -1.7857142857146e-02; +Q11_15 = -1.7857142857146e-02; +Q11_16 = 2.5974025974031e-03; +Q11_16 = 2.5974025974031e-03; +Q11_17 = -1.8037518037522e-04; +Q11_17 = -1.8037518037522e-04; +for i = 1:BP + for j = 1:BP + Q(i,j) = eval(['Q' num2str(i-1) '_' num2str(j-1)]); + Q(N+1-i,N+1-j) = -eval(['Q' num2str(i-1) '_' num2str(j-1)]); + end +end +%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%% Difference operator %% +D1 = H\Q; +%%%%%%%%%%%%%%%%%%%%%%%%%%% \ No newline at end of file