Mercurial > repos > public > sbplib
diff +sbp/+implementations/d2_noneq_variable_12.m @ 1325:1b0f2415237f feature/D2_boundary_opt
Add variable coefficient boundary-optimized second derivatives.
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Sun, 13 Feb 2022 19:32:34 +0100 |
| parents | |
| children | c2d716c4f1ed |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d2_noneq_variable_12.m Sun Feb 13 19:32:34 2022 +0100 @@ -0,0 +1,401 @@ +function [H, HI, D1, D2, DI] = d2_noneq_variable_12(N, h, options) + % N: Number of grid points + % h: grid spacing + % options: struct containing options for constructing the operator + % current options are: + % options.stencil_type ('minimal','nonminimal','wide') + % options.AD ('upwind', 'op') + + % BP: Number of boundary points + % order: Accuracy of interior stencil + BP = 12; + order = 12; + + %%%% Norm matrix %%%%%%%% + P = zeros(BP, 1); + 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 + + Hv = ones(N, 1); + Hv(1:BP) = P; + Hv(end - BP + 1:end) = flip(P); + Hv = h * Hv; + H = spdiags(Hv, 0, N, N); + HI = spdiags(1 ./ Hv, 0, N, N); + %%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Q matrix %%%%%%%%%%% + + % interior stencil + 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]; + 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_1 = 7.8128092862319e-04; + Q10_2 = -1.6271304373071e-03; + Q10_3 = 3.2632650250470e-03; + Q10_4 = -6.1239492854797e-03; + Q10_5 = 8.4387621360184e-03; + Q10_6 = 4.3175367549700e-03; + Q10_7 = -6.9373143801571e-02; + Q10_8 = 2.6369594097148e-01; + Q10_9 = -8.5631774953989e-01; + Q10_10 = 0.0000000000000e+00; + Q10_11 = 8.5712580212095e-01; + Q10_12 = -2.6785714285718e-01; + Q10_13 = 7.9365079365093e-02; + Q10_14 = -1.7857142857146e-02; + Q10_15 = 2.5974025974031e-03; + Q10_16 = -1.8037518037522e-04; + Q10_17 = 0.0000000000000e+00; + Q11_0 = -3.1501105449828e-05; + Q11_1 = 6.0444181254875e-05; + Q11_2 = -1.7066984372933e-05; + Q11_3 = -1.3097937809499e-04; + Q11_4 = 4.5892226603067e-04; + Q11_5 = -8.5964292632428e-04; + Q11_6 = -1.1092605832824e-03; + Q11_7 = 1.6606121869177e-02; + Q11_8 = -7.8759594625702e-02; + Q11_9 = 2.6769768119702e-01; + Q11_10 = -8.5712580212095e-01; + Q11_11 = 0.0000000000000e+00; + Q11_12 = 8.5714285714289e-01; + Q11_13 = -2.6785714285718e-01; + Q11_14 = 7.9365079365093e-02; + Q11_15 = -1.7857142857146e-02; + Q11_16 = 2.5974025974031e-03; + 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 + + %%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%% Undivided difference operators %%%% + % Closed with zeros at the first boundary nodes. + m = N; + + DD_6 = (diag(ones(m - 3, 1), 3) - 6 * diag(ones(m - 2, 1), 2) + 15 * diag(ones(m - 1, 1), 1) - 20 * diag(ones(m, 1), 0) + 15 * diag(ones(m - 1, 1), -1) - 6 * diag(ones(m - 2, 1), -2) + diag(ones(m - 3, 1), -3)); + DD_6(1:9, 1:12) = [0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0.70504538217624485042e1 -0.15481773512478815978e2 0.15439628908176388584e2 -0.11355022328644135767e2 0.62553652361133177256e1 -0.23565505827757565576e1 0.44789845784655348861e0 0 0 0 0 0; 0 0.84178325750049836493e0 -0.30776567832262484696e1 0.55480476621608890820e1 -0.66451562589578036671e1 0.54681670852508283306e1 -0.26873907470793164527e1 0.55220578435115281189e0 0 0 0 0; 0 0 0.36124410255434790612e0 -0.18841869963752696615e1 0.49068751071627186058e1 -0.79710602635488243814e1 0.75771246506939812399e1 -0.36661050678180486359e1 0.67610846733109492706e0 0 0 0; 0 0 0 0.31883568286286899671e0 -0.22216567686182446798e1 0.72341822309820977868e1 -0.12215760254444741754e2 0.10698736280837039597e2 -0.46222617037529123987e1 0.80792453213389245238e0 0 0; 0 0 0 0 0.45451605753051033631e0 -0.36739526096585069959e1 0.11306936594922967126e2 -0.16769946247690595362e2 0.13179014442979029716e2 -0.54150410306154005497e1 0.91847279253199572926e0 0; 0 0 0 0 0 0.77355004999207313207e0 -0.54143198515959566716e1 0.14230335022999000858e2 -0.19303948318184081752e2 0.14605034473743418451e2 -0.58729419174319386168e1 0.98229054047748459948e0; ]; + DD_6(m - 8:m, m - 11:m) = [0.98229054047748459948e0 -0.58729419174319386168e1 0.14605034473743418451e2 -0.19303948318184081752e2 0.14230335022999000858e2 -0.54143198515959566716e1 0.77355004999207313207e0 0 0 0 0 0; 0 0.91847279253199572926e0 -0.54150410306154005497e1 0.13179014442979029716e2 -0.16769946247690595362e2 0.11306936594922967126e2 -0.36739526096585069959e1 0.45451605753051033631e0 0 0 0 0; 0 0 0.80792453213389245238e0 -0.46222617037529123987e1 0.10698736280837039597e2 -0.12215760254444741754e2 0.72341822309820977868e1 -0.22216567686182446798e1 0.31883568286286899671e0 0 0 0; 0 0 0 0.67610846733109492706e0 -0.36661050678180486359e1 0.75771246506939812399e1 -0.79710602635488243814e1 0.49068751071627186058e1 -0.18841869963752696615e1 0.36124410255434790612e0 0 0; 0 0 0 0 0.55220578435115281189e0 -0.26873907470793164527e1 0.54681670852508283306e1 -0.66451562589578036671e1 0.55480476621608890820e1 -0.30776567832262484696e1 0.84178325750049836493e0 0; 0 0 0 0 0 0.44789845784655348861e0 -0.23565505827757565576e1 0.62553652361133177256e1 -0.11355022328644135767e2 0.15439628908176388584e2 -0.15481773512478815978e2 0.70504538217624485042e1; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_6 = sparse(DD_6); + + DD_7 = (-diag(ones(m - 4, 1), -4) + 7 * diag(ones(m - 3, 1), -3) - 21 * diag(ones(m - 2, 1), -2) + 35 * diag(ones(m - 1, 1), -1) - 35 * diag(ones(m, 1), 0) + 21 * diag(ones(m - 1, 1), 1) - 7 * diag(ones(m - 2, 1), 2) + diag(ones(m - 3, 1), 3)); + DD_7(1:10, 1:13) = [0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; -0.70504538217624585762e1 0.16323556769979337662e2 -0.18517285691402663507e2 0.16903069990805048996e2 -0.12900521495071139822e2 0.78247176680265960663e1 -0.31352892049258744203e1 0.55220578435115360075e0 0 0 0 0 0; 0 -0.77136636008199086135e0 0.31512297676528003033e1 -0.68105129732006589582e1 0.10585680229488408490e2 -0.12315008394369015993e2 0.94058676147776075845e1 -0.38654404904580696832e1 0.61955060619091911792e0 0 0 0 0; 0 0 -0.32268062071305431283e0 0.19678458985349116625e1 -0.63675477999083230026e1 0.13582054493165302513e2 -0.17679957518285956226e2 0.12831367737363170226e2 -0.47327592713176644894e1 0.72167708116161363042e0 0 0 0; 0 0 0 -0.28975994984220671445e0 0.24321233335964448997e1 -0.99133841479577475282e1 0.21377580445278298070e2 -0.24963717988619759059e2 0.16177915963135193396e2 -0.56554717249372471666e1 0.83471406934702410357e0 0 0; 0 0 0 0 -0.43028213606032104517e0 0.42103703770684731501e1 -0.15829711232892153976e2 0.29347405933458541883e2 -0.30751033700284402671e2 0.18952643607153901924e2 -0.64293095477239701048e1 0.92991669927993083983e0 0; 0 0 0 0 0 -0.76177813656089336989e0 0.63167064935286161169e1 -0.19922469032198601202e2 0.33781909556822143066e2 -0.34078413772067976385e2 0.20555296711011785159e2 -0.68760337833423921963e1 0.98478196280731881078e0; ]; + DD_7(m - 8:m, m - 12:m) = [-0.98478196280731881078e0 0.68760337833423921963e1 -0.20555296711011785159e2 0.34078413772067976385e2 -0.33781909556822143066e2 0.19922469032198601202e2 -0.63167064935286161169e1 0.76177813656089336989e0 0 0 0 0 0; 0 -0.92991669927993083983e0 0.64293095477239701048e1 -0.18952643607153901924e2 0.30751033700284402671e2 -0.29347405933458541883e2 0.15829711232892153976e2 -0.42103703770684731501e1 0.43028213606032104517e0 0 0 0 0; 0 0 -0.83471406934702410357e0 0.56554717249372471666e1 -0.16177915963135193396e2 0.24963717988619759059e2 -0.21377580445278298070e2 0.99133841479577475282e1 -0.24321233335964448997e1 0.28975994984220671445e0 0 0 0; 0 0 0 -0.72167708116161363042e0 0.47327592713176644894e1 -0.12831367737363170226e2 0.17679957518285956226e2 -0.13582054493165302513e2 0.63675477999083230026e1 -0.19678458985349116625e1 0.32268062071305431283e0 0 0; 0 0 0 0 -0.61955060619091911792e0 0.38654404904580696832e1 -0.94058676147776075845e1 0.12315008394369015993e2 -0.10585680229488408490e2 0.68105129732006589582e1 -0.31512297676528003033e1 0.77136636008199086135e0 0; 0 0 0 0 0 -0.55220578435115360075e0 0.31352892049258744203e1 -0.78247176680265960663e1 0.12900521495071139822e2 -0.16903069990805048996e2 0.18517285691402663507e2 -0.16323556769979337662e2 0.70504538217624585762e1; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_7 = sparse(DD_7); + + DD_8 = (diag(ones(m - 4, 1), 4) - 8 * diag(ones(m - 3, 1), 3) + 28 * diag(ones(m - 2, 1), 2) - 56 * diag(ones(m - 1, 1), 1) + 70 * diag(ones(m, 1), 0) - 56 * diag(ones(m - 1, 1), -1) + 28 * diag(ones(m - 2, 1), -2) - 8 * diag(ones(m - 3, 1), -3) + diag(ones(m - 4, 1), -4)); + DD_8(1:10, 1:14) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0.70504538217624673892e1 -0.17094923130061349891e2 0.21668515459055490895e2 -0.23713582964005737596e2 0.23486201724559577669e2 -0.20139726062395637234e2 0.12541156819703497681e2 -0.44176462748092288060e1 0.61955060619091989235e0 0 0 0 0 0; 0 0.71430915910501867497e0 -0.32169486987428931518e1 0.81290324137612243914e1 -0.15699214646211457760e2 0.23981482952559874949e2 -0.25082313639406953559e2 0.15461761961832278733e2 -0.49564048495273529434e1 0.66829534663026066516e0 0 0 0 0; 0 0 0.29213206789956748018e0 -0.20438756549159061384e1 0.79665959467484077329e1 -0.21271097908734749058e2 0.35359915036571912453e2 -0.34216980632968453935e2 0.18931037085270657958e2 -0.57734166492929090433e1 0.75569070942147255138e0 0 0 0; 0 0 0 0.26637215914130374043e0 -0.26313682263734604148e1 0.12983764630211125301e2 -0.34204128712445276912e2 0.49927435977239518119e2 -0.43141109235027182388e2 0.22621886899748988667e2 -0.66777125547761928286e1 0.85485906228117671607e0 0 0; 0 0 0 0 0.41007336094698233166e0 -0.47386249189431794318e1 0.21106281643856205301e2 -0.46955849493533667013e2 0.61502067400568805342e2 -0.50540382952410405130e2 0.25717238190895880419e2 -0.74393335942394467187e1 0.93853036285882489957e0 0; 0 0 0 0 0 0.75161513192516571838e0 -0.72190931354612755621e1 0.26563292042931468269e2 -0.54051055290915428906e2 0.68156827544135952769e2 -0.54814124562698093757e2 0.27504135133369568785e2 -0.78782557024585504862e1 0.98665883917119316896e0; ]; + DD_8(m - 9:m, m - 13:m) = [0.98665883917119316896e0 -0.78782557024585504862e1 0.27504135133369568785e2 -0.54814124562698093757e2 0.68156827544135952769e2 -0.54051055290915428906e2 0.26563292042931468269e2 -0.72190931354612755621e1 0.75161513192516571838e0 0 0 0 0 0; 0 0.93853036285882489957e0 -0.74393335942394467187e1 0.25717238190895880419e2 -0.50540382952410405130e2 0.61502067400568805342e2 -0.46955849493533667013e2 0.21106281643856205301e2 -0.47386249189431794318e1 0.41007336094698233166e0 0 0 0 0; 0 0 0.85485906228117671607e0 -0.66777125547761928286e1 0.22621886899748988667e2 -0.43141109235027182388e2 0.49927435977239518119e2 -0.34204128712445276912e2 0.12983764630211125301e2 -0.26313682263734604148e1 0.26637215914130374043e0 0 0 0; 0 0 0 0.75569070942147255138e0 -0.57734166492929090433e1 0.18931037085270657958e2 -0.34216980632968453935e2 0.35359915036571912453e2 -0.21271097908734749058e2 0.79665959467484077329e1 -0.20438756549159061384e1 0.29213206789956748018e0 0 0; 0 0 0 0 0.66829534663026066516e0 -0.49564048495273529434e1 0.15461761961832278733e2 -0.25082313639406953559e2 0.23981482952559874949e2 -0.15699214646211457760e2 0.81290324137612243914e1 -0.32169486987428931518e1 0.71430915910501867497e0 0; 0 0 0 0 0 0.61955060619091989235e0 -0.44176462748092288060e1 0.12541156819703497681e2 -0.20139726062395637234e2 0.23486201724559577669e2 -0.23713582964005737596e2 0.21668515459055490895e2 -0.17094923130061349891e2 0.70504538217624673892e1; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_8 = sparse(DD_8); + + DD_9 = (-diag(ones(m - 5, 1), -5) + 9 * diag(ones(m - 4, 1), -4) - 36 * diag(ones(m - 3, 1), -3) + 84 * diag(ones(m - 2, 1), -2) - 126 * diag(ones(m - 1, 1), -1) + 126 * diag(ones(m, 1), 0) - 84 * diag(ones(m - 1, 1), 1) + 36 * diag(ones(m - 2, 1), 2) - 9 * diag(ones(m - 3, 1), 3) + diag(ones(m - 4, 1), 4)); + DD_9(1:11, 1:15) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; -0.70504538217624752230e1 0.17809232289166388354e2 -0.24885464157798411697e2 0.31842615377766997368e2 -0.39185416370771078968e2 0.44121209014955561206e2 -0.37623470459110493043e2 0.19879408236641529627e2 -0.55759554557182790312e1 0.66829534663026140771e0 0 0 0 0 0; 0 -0.66695396914459764542e0 0.32764459415493351213e1 -0.94984942131335544919e1 0.22096883550779531311e2 -0.42252556942280502185e2 0.56435205688665645507e2 -0.46385285885496836199e2 0.22303821822873088245e2 -0.60146581196723459865e1 0.70559212586023632254e0 0 0 0 0; 0 0 -0.26728718140337933088e0 0.21137687176984662199e1 -0.96966416347745772182e1 0.31341597391980823551e2 -0.63647847065829442415e2 0.76988206424179021354e2 -0.56793111255811973873e2 0.25980374921818090695e2 -0.68012163847932529624e1 0.78215606693622397877e0 0 0 0; 0 0 0 -0.24708804973573377055e0 0.28212553166921198174e1 -0.16439370707786854975e2 0.51306193068667915368e2 -0.89869384759031132614e2 0.97067495778811160373e2 -0.67865660699246966000e2 0.30049706496492867728e2 -0.76937315605305904447e1 0.87058511566721451662e0 0 0; 0 0 0 0 -0.39286387086790423439e0 0.52598319320161875843e1 -0.27136647827815121102e2 0.70433774240300500520e2 -0.11070372132102384962e3 0.11371586164292341154e3 -0.77151714572687641258e2 0.33477001174077510234e2 -0.84467732657294240961e1 0.94525186880633042479e0 0; 0 0 0 0 0 -0.74268863896636254047e0 0.81214797773939350074e1 -0.34152804055197602060e2 0.81076582936373143359e2 -0.12268228957944471498e3 0.12333178026607071095e3 -0.82512405400108706356e2 0.35452150661063477188e2 -0.88799295525407385207e1 0.98812358535685795524e0; ]; + DD_9(m - 9:m, m - 14:m) = [-0.98812358535685795524e0 0.88799295525407385207e1 -0.35452150661063477188e2 0.82512405400108706356e2 -0.12333178026607071095e3 0.12268228957944471498e3 -0.81076582936373143359e2 0.34152804055197602060e2 -0.81214797773939350074e1 0.74268863896636254047e0 0 0 0 0 0; 0 -0.94525186880633042479e0 0.84467732657294240961e1 -0.33477001174077510234e2 0.77151714572687641258e2 -0.11371586164292341154e3 0.11070372132102384962e3 -0.70433774240300500520e2 0.27136647827815121102e2 -0.52598319320161875843e1 0.39286387086790423439e0 0 0 0 0; 0 0 -0.87058511566721451662e0 0.76937315605305904447e1 -0.30049706496492867728e2 0.67865660699246966000e2 -0.97067495778811160373e2 0.89869384759031132614e2 -0.51306193068667915368e2 0.16439370707786854975e2 -0.28212553166921198174e1 0.24708804973573377055e0 0 0 0; 0 0 0 -0.78215606693622397877e0 0.68012163847932529624e1 -0.25980374921818090695e2 0.56793111255811973873e2 -0.76988206424179021354e2 0.63647847065829442415e2 -0.31341597391980823551e2 0.96966416347745772182e1 -0.21137687176984662199e1 0.26728718140337933088e0 0 0; 0 0 0 0 -0.70559212586023632254e0 0.60146581196723459865e1 -0.22303821822873088245e2 0.46385285885496836199e2 -0.56435205688665645507e2 0.42252556942280502185e2 -0.22096883550779531311e2 0.94984942131335544919e1 -0.32764459415493351213e1 0.66695396914459764542e0 0; 0 0 0 0 0 -0.66829534663026140771e0 0.55759554557182790312e1 -0.19879408236641529627e2 0.37623470459110493043e2 -0.44121209014955561206e2 0.39185416370771078968e2 -0.31842615377766997368e2 0.24885464157798411697e2 -0.17809232289166388354e2 0.70504538217624752230e1; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_9 = sparse(DD_9); + + DD_10 = (diag(ones(m - 5, 1), 5) - 10 * diag(ones(m - 4, 1), 4) + 45 * diag(ones(m - 3, 1), 3) - 120 * diag(ones(m - 2, 1), 2) + 210 * diag(ones(m - 1, 1), 1) - 252 * diag(ones(m, 1), 0) + 210 * diag(ones(m - 1, 1), -1) - 120 * diag(ones(m - 2, 1), -2) + 45 * diag(ones(m - 3, 1), -3) - 10 * diag(ones(m - 4, 1), -4) + diag(ones(m - 5, 1), -5)); + DD_10(1:11, 1:16) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0.70504538217624822734e1 -0.18476186258311004476e2 0.28161910099347774980e2 -0.41341109590900593201e2 0.61282299921550671561e2 -0.86373765957236149764e2 0.94058676147776232608e2 -0.66264694122138432090e2 0.27879777278591395156e2 -0.66829534663026140771e1 0.70559212586023702812e0 0 0 0 0 0; 0 0.62689419647750179604e0 -0.33308831364980034247e1 0.10914786591113557343e2 -0.29883886064802801254e2 0.69173811658980680918e2 -0.11287041137733129101e3 0.11596321471374209050e3 -0.74346072742910294151e2 0.30073290598361729932e2 -0.70559212586023632254e1 0.73517682146919258207e0 0 0 0 0; 0 0 0.24665297575614270036e0 -0.21786018467637256442e1 0.11551532389532584562e2 -0.44092342567416599454e2 0.10607974510971573736e3 -0.15397641284835804271e3 0.14198277813952993468e3 -0.86601249739393635650e2 0.34006081923966264812e2 -0.78215606693622397877e1 0.80337713279357912969e0 0 0 0; 0 0 0 0.23087142251463535911e0 -0.30031734426541638399e1 0.20275063024062368195e2 -0.73294561526668450525e2 0.14978230793171855436e3 -0.19413499155762232075e3 0.16966415174811741500e3 -0.10016568832164289243e3 0.38468657802652952223e2 -0.87058511566721451662e1 0.88321407619404757308e0 0 0; 0 0 0 0 0.37796280007363075810e0 -0.57748488744870271355e1 0.33920809784768901377e2 -0.10061967748614357217e3 0.18450620220170641603e3 -0.22743172328584682309e3 0.19287928643171910314e3 -0.11159000391359170078e3 0.42233866328647120480e2 -0.94525186880633042479e1 0.95064470121725563376e0 0; 0 0 0 0 0 0.73474076934244039806e0 -0.90238664193265944527e1 0.42691005068997002575e2 -0.11582368990910449051e3 0.20447048263240785831e3 -0.24666356053214142190e3 0.20628101350027176589e3 -0.11817383553687825729e3 0.44399647762703692603e2 -0.98812358535685795524e1 0.98929851729658394154e0; ]; + DD_10(m - 10:m, m - 15:m) = [0.98929851729658394154e0 -0.98812358535685795524e1 0.44399647762703692603e2 -0.11817383553687825729e3 0.20628101350027176589e3 -0.24666356053214142190e3 0.20447048263240785831e3 -0.11582368990910449051e3 0.42691005068997002575e2 -0.90238664193265944527e1 0.73474076934244039806e0 0 0 0 0 0; 0 0.95064470121725563376e0 -0.94525186880633042479e1 0.42233866328647120480e2 -0.11159000391359170078e3 0.19287928643171910314e3 -0.22743172328584682309e3 0.18450620220170641603e3 -0.10061967748614357217e3 0.33920809784768901377e2 -0.57748488744870271355e1 0.37796280007363075810e0 0 0 0 0; 0 0 0.88321407619404757308e0 -0.87058511566721451662e1 0.38468657802652952223e2 -0.10016568832164289243e3 0.16966415174811741500e3 -0.19413499155762232075e3 0.14978230793171855436e3 -0.73294561526668450525e2 0.20275063024062368195e2 -0.30031734426541638399e1 0.23087142251463535911e0 0 0 0; 0 0 0 0.80337713279357912969e0 -0.78215606693622397877e1 0.34006081923966264812e2 -0.86601249739393635650e2 0.14198277813952993468e3 -0.15397641284835804271e3 0.10607974510971573736e3 -0.44092342567416599454e2 0.11551532389532584562e2 -0.21786018467637256442e1 0.24665297575614270036e0 0 0; 0 0 0 0 0.73517682146919258207e0 -0.70559212586023632254e1 0.30073290598361729932e2 -0.74346072742910294151e2 0.11596321471374209050e3 -0.11287041137733129101e3 0.69173811658980680918e2 -0.29883886064802801254e2 0.10914786591113557343e2 -0.33308831364980034247e1 0.62689419647750179604e0 0; 0 0 0 0 0 0.70559212586023702812e0 -0.66829534663026140771e1 0.27879777278591395156e2 -0.66264694122138432090e2 0.94058676147776232608e2 -0.86373765957236149764e2 0.61282299921550671561e2 -0.41341109590900593201e2 0.28161910099347774980e2 -0.18476186258311004476e2 0.70504538217624822734e1; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_10 = sparse(DD_10); + + %[-1, 11, -55, 165, -330, 462, -462, 330, -165, 55, -11, 1, 0] + DD_11 = (-diag(ones(m - 6, 1), -6) + 11 * diag(ones(m - 5, 1), -5) - 55 * diag(ones(m - 4, 1), -4) + 165 * diag(ones(m - 3, 1), -3) - 330 * diag(ones(m - 2, 1), -2) + 462 * diag(ones(m - 1, 1), -1) - 462 * diag(ones(m, 1), 0) + 330 * diag(ones(m - 1, 1), 1) - 165 * diag(ones(m - 2, 1), 2) + 55 * diag(ones(m - 3, 1), 3) - 11 * diag(ones(m - 4, 1), 4) + diag(ones(m - 5, 1), 5)); + DD_11(1:12, 1:17) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; -0.70504538217624886828e1 0.19103080454788523638e2 -0.31492793235845807035e2 0.52255896182014198048e2 -0.91166185986353555693e2 0.15554757761621697209e3 -0.20692908752510771174e3 0.18222790883588068825e3 -0.10222585002150178224e3 0.36756244064664377424e2 -0.77615133844626073094e1 0.73517682146919325041e0 0 0 0 0 0; 0 -0.59247568548574727441e0 0.33811178542850964393e1 -0.12374519230919224800e2 0.39160480492663455349e2 -0.10704747746062038309e3 0.20692908752510736686e3 -0.25511907237023259909e3 0.20445170004300330891e3 -0.11026873219399300975e3 0.38807566922312997740e2 -0.80869450361611184028e1 0.75926914003985711330e0 0 0 0 0; 0 0 -0.22922038452398770506e0 0.22391799149842653766e1 -0.13526028855965613027e2 0.59818136859120542509e2 -0.16669674231526758728e3 0.28229009022198974497e3 -0.31236211190696585630e3 0.23815343678333249804e3 -0.12468896705454297098e3 0.43018583681492318832e2 -0.88371484607293704266e1 0.82079151707601599340e0 0 0 0; 0 0 0 -0.21701391114437616116e0 0.31781915325611402199e1 -0.24486327517266714070e2 0.10078002209916911947e3 -0.23537219817841487113e3 0.35591415118897425470e3 -0.37326113384585831300e3 0.27545564288451795418e3 -0.14105174527639415815e3 0.47882181361696798414e2 -0.97153548381345233039e1 0.89358450029368883150e0 0 0; 0 0 0 0 -0.36488508570871329747e0 0.62843543720560487741e1 -0.41458767514717546128e2 0.13835205654344741174e3 -0.28993831774553865375e3 0.41695815935738584232e3 -0.42433443014978202692e3 0.30687251076237717714e3 -0.15485750987170610843e3 0.51988852784348173364e2 -0.10457091713389811971e2 0.95506826122820715686e0 0; 0 0 0 0 0 -0.72758579432539352184e0 0.99262530612592538979e1 -0.52177895084329669814e2 0.15925757362501867446e3 -0.32131075842235520591e3 0.45221652764225927349e3 -0.45381822970059788496e3 0.32497804772641520756e3 -0.16279870846324687288e3 0.54346797194627187538e2 -0.10882283690262423357e2 0.99026190553785350209e0; ]; + DD_11(m - 10:m, m - 16:m) = [-0.99026190553785350209e0 0.10882283690262423357e2 -0.54346797194627187538e2 0.16279870846324687288e3 -0.32497804772641520756e3 0.45381822970059788496e3 -0.45221652764225927349e3 0.32131075842235520591e3 -0.15925757362501867446e3 0.52177895084329669814e2 -0.99262530612592538979e1 0.72758579432539352184e0 0 0 0 0 0; 0 -0.95506826122820715686e0 0.10457091713389811971e2 -0.51988852784348173364e2 0.15485750987170610843e3 -0.30687251076237717714e3 0.42433443014978202692e3 -0.41695815935738584232e3 0.28993831774553865375e3 -0.13835205654344741174e3 0.41458767514717546128e2 -0.62843543720560487741e1 0.36488508570871329747e0 0 0 0 0; 0 0 -0.89358450029368883150e0 0.97153548381345233039e1 -0.47882181361696798414e2 0.14105174527639415815e3 -0.27545564288451795418e3 0.37326113384585831300e3 -0.35591415118897425470e3 0.23537219817841487113e3 -0.10078002209916911947e3 0.24486327517266714070e2 -0.31781915325611402199e1 0.21701391114437616116e0 0 0 0; 0 0 0 -0.82079151707601599340e0 0.88371484607293704266e1 -0.43018583681492318832e2 0.12468896705454297098e3 -0.23815343678333249804e3 0.31236211190696585630e3 -0.28229009022198974497e3 0.16669674231526758728e3 -0.59818136859120542509e2 0.13526028855965613027e2 -0.22391799149842653766e1 0.22922038452398770506e0 0 0; 0 0 0 0 -0.75926914003985711330e0 0.80869450361611184028e1 -0.38807566922312997740e2 0.11026873219399300975e3 -0.20445170004300330891e3 0.25511907237023259909e3 -0.20692908752510736686e3 0.10704747746062038309e3 -0.39160480492663455349e2 0.12374519230919224800e2 -0.33811178542850964393e1 0.59247568548574727441e0 0; 0 0 0 0 0 -0.73517682146919325041e0 0.77615133844626073094e1 -0.36756244064664377424e2 0.10222585002150178224e3 -0.18222790883588068825e3 0.20692908752510771174e3 -0.15554757761621697209e3 0.91166185986353555693e2 -0.52255896182014198048e2 0.31492793235845807035e2 -0.19103080454788523638e2 0.70504538217624886828e1; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_11 = sparse(DD_11); + + %[1, -12, 66, -220, 495, -792, 924, -792, 495, -220, 66,-12, 1] + DD_12 = (diag(ones(m - 6, 1), 6) - 12 * diag(ones(m - 5, 1), 5) + 66 * diag(ones(m - 4, 1), 4) - 220 * diag(ones(m - 3, 1), 3) + 495 * diag(ones(m - 2, 1), 2) - 792 * diag(ones(m - 1, 1), 1) + 924 * diag(ones(m, 1), 0) - 792 * diag(ones(m - 1, 1), -1) + 495 * diag(ones(m - 2, 1), -2) - 220 * diag(ones(m - 3, 1), -3) + 66 * diag(ones(m - 4, 1), -4) - 12 * diag(ones(m - 5, 1), -5) + diag(ones(m - 6, 1), -6)); + DD_12(1:12, 1:18) = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0.70504538217624945581e1 -0.19695556140274287325e2 0.34873911090130932535e2 -0.64630415412933476707e2 0.13032666647901711965e3 -0.26259505507683757401e3 0.41385817505021542347e3 -0.43734698120611365179e3 0.30667755006450534672e3 -0.14702497625865750970e3 0.46569080306775643856e2 -0.88221218576303190049e1 0.75926914003985774602e0 0 0 0 0 0; 0 0.56252054413792412864e0 -0.34278021535209417538e1 0.13874841106233553089e2 -0.50022717612825328606e2 0.15842900977125025845e3 -0.35473557861446977176e3 0.51023814474046519819e3 -0.49068408010320794140e3 0.33080619658197902926e3 -0.15523026768925199096e3 0.48521670216966710417e2 -0.91112296804782853596e1 0.77929289272158630803e0 0 0 0 0; 0 0 0.21428192906429919385e0 -0.22961219278627835050e1 0.15615594467315483385e2 -0.78810282483468544556e2 0.25004511347290138092e3 -0.48392586895198241994e3 0.62472422381393171260e3 -0.57156824827999799529e3 0.37406690116362891293e3 -0.17207433472596927533e3 0.53022890764376222560e2 -0.98494982049121919208e1 0.83534896297519895228e0 0 0 0; 0 0 0 0.20501361895561890151e0 -0.33471539032596363119e1 0.29069145008244604383e2 -0.13437336279889215930e3 0.35305829726762230670e3 -0.61013854489538443663e3 0.74652226769171662600e3 -0.66109354292284309003e3 0.42315523582918247446e3 -0.19152872544678719366e3 0.58292129028807139823e2 -0.10723014003524265978e2 0.90225552616201164306e0 0 0; 0 0 0 0 0.35327850260839005288e0 -0.67888982569607603271e1 0.49750521017661055353e2 -0.18446940872459654898e3 0.43490747661830798063e3 -0.71478541604123287255e3 0.84866886029956405384e3 -0.73649402582970522515e3 0.46457252961511832529e3 -0.20795541113739269345e3 0.62742550280338871828e2 -0.11460819134738485882e2 0.95876279102790935799e0 0; 0 0 0 0 0 0.72108566182178027310e0 -0.10828639703191913343e2 0.62613474101195603777e2 -0.21234343150002489927e3 0.48196613763353280887e3 -0.77522833310101589741e3 0.90763645940119576992e3 -0.77994731454339649814e3 0.48839612538974061864e3 -0.21738718877850875015e3 0.65293702141574540142e2 -0.11883142866454242025e2 0.99106616353107873319e0; ]; + DD_12(m - 11:m, m - 17:m) = [0.99106616353107873319e0 -0.11883142866454242025e2 0.65293702141574540142e2 -0.21738718877850875015e3 0.48839612538974061864e3 -0.77994731454339649814e3 0.90763645940119576992e3 -0.77522833310101589741e3 0.48196613763353280887e3 -0.21234343150002489927e3 0.62613474101195603777e2 -0.10828639703191913343e2 0.72108566182178027310e0 0 0 0 0 0; 0 0.95876279102790935799e0 -0.11460819134738485882e2 0.62742550280338871828e2 -0.20795541113739269345e3 0.46457252961511832529e3 -0.73649402582970522515e3 0.84866886029956405384e3 -0.71478541604123287255e3 0.43490747661830798063e3 -0.18446940872459654898e3 0.49750521017661055353e2 -0.67888982569607603271e1 0.35327850260839005288e0 0 0 0 0; 0 0 0.90225552616201164306e0 -0.10723014003524265978e2 0.58292129028807139823e2 -0.19152872544678719366e3 0.42315523582918247446e3 -0.66109354292284309003e3 0.74652226769171662600e3 -0.61013854489538443663e3 0.35305829726762230670e3 -0.13437336279889215930e3 0.29069145008244604383e2 -0.33471539032596363119e1 0.20501361895561890151e0 0 0 0; 0 0 0 0.83534896297519895228e0 -0.98494982049121919208e1 0.53022890764376222560e2 -0.17207433472596927533e3 0.37406690116362891293e3 -0.57156824827999799529e3 0.62472422381393171260e3 -0.48392586895198241994e3 0.25004511347290138092e3 -0.78810282483468544556e2 0.15615594467315483385e2 -0.22961219278627835050e1 0.21428192906429919385e0 0 0; 0 0 0 0 0.77929289272158630803e0 -0.91112296804782853596e1 0.48521670216966710417e2 -0.15523026768925199096e3 0.33080619658197902926e3 -0.49068408010320794140e3 0.51023814474046519819e3 -0.35473557861446977176e3 0.15842900977125025845e3 -0.50022717612825328606e2 0.13874841106233553089e2 -0.34278021535209417538e1 0.56252054413792412864e0 0; 0 0 0 0 0 0.75926914003985774602e0 -0.88221218576303190049e1 0.46569080306775643856e2 -0.14702497625865750970e3 0.30667755006450534672e3 -0.43734698120611365179e3 0.41385817505021542347e3 -0.26259505507683757401e3 0.13032666647901711965e3 -0.64630415412933476707e2 0.34873911090130932535e2 -0.19695556140274287325e2 0.70504538217624945581e1; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; ]; + DD_12 = sparse(DD_12); + %%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Difference operators %% + D1 = H \ Q; + + % Helper functions for constructing D2(c) + % TODO: Consider changing sparse(diag(...)) to spdiags(....) + + % Minimal 13 point stencil width + function D2 = D2_fun_minimal(c) + % Here we add variable diffusion + C1 = sparse(diag(c)); + C2 = 1/2 * diag(ones(m - 1, 1), -1) + 1/2 * diag(ones(m, 1), 0); C2(1, 2) = 1/2; + C3 = 5/18 * diag(ones(m - 1, 1), -1) + 5/18 * diag(ones(m - 1, 1), 1) + 4/9 * diag(ones(m, 1), 0); C3(1, 3) = 5/18; C3(m, m - 2) = 5/18; + C4 = 5/28 * diag(ones(m - 2, 1), -2) + 9/28 * diag(ones(m - 1, 1), -1) + 5/28 * diag(ones(m - 1, 1), 1) + 9/28 * diag(ones(m, 1), 0); C4(2, 4) = 5/28; C4(1, 3) = 5/28; C4(1, 4) = 9/28; C4(m, m - 2) = 5/28; C4(m, m - 3) = 5/28; + C5 = 1/7 * diag(ones(m - 2, 1), -2) + 1/7 * diag(ones(m - 2, 1), 2) + 8/35 * diag(ones(m - 1, 1), -1) + 8/35 * diag(ones(m - 1, 1), 1) + 9/35 * diag(ones(m, 1), 0); C5(2, 4) = 1/7; C5(1, 3) = 1/7; C5(1, 4) = 1/7; C5(1, 5) = 8/35; C5(2, 5) = 1/7; C5(m - 1, m - 4) = 1/7; C5(m, m - 4) = 8/35; C5(m, m - 3) = 1/7; + C6 = 1/6 * diag(ones(m - 3, 1), -3) + 1/6 * diag(ones(m - 2, 1), -2) + 1/6 * diag(ones(m - 2, 1), 2) + 1/6 * diag(ones(m - 1, 1), -1) + 1/6 * diag(ones(m - 1, 1), 1) + 1/6 * diag(ones(m, 1), 0); + C6(1, 4) = 1/6; C6(1, 5) = 1/6; C6(1, 6) = 1/6; C6(2, 5) = 1/6; C6(2, 6) = 1/6; C6(3, 6) = 1/6; C6(m - 1, m - 5) = 1/6; C6(m, m - 4) = 1/6; C6(m, m - 5) = 1/6; + + C2 = sparse(diag(C2 * c)); + C3 = sparse(diag(C3 * c)); + C4 = sparse(diag(C4 * c)); + C5 = sparse(diag(C5 * c)); + C6 = sparse(diag(C6 * c)); + + % Remainder term added to wide second derivative operator + R = (1/30735936 / h) * transpose(DD_12) * C1 * DD_12 + (1/6403320 / h) * transpose(DD_11) * C2 * DD_11 + (1/1293600 / h) * transpose(DD_10) * C3 * DD_10 + (1/249480 / h) * transpose(DD_9) * C4 * DD_9 + (1/44352 / h) * transpose(DD_8) * C5 * DD_8 + (1/6468 / h) * transpose(DD_7) * C6 * DD_7; + D2 = D1 * C1 * D1 - H \ R; + end + + % Few additional grid point in interior stencil cmp. to minimal + function D2 = D2_fun_nonminimal(c) + % Here we add variable diffusion + C1 = sparse(diag(c)); + C2 = 1/2 * diag(ones(m - 1, 1), -1) + 1/2 * diag(ones(m, 1), 0); C2(1, 2) = 1/2; + C3 = 5/18 * diag(ones(m - 1, 1), -1) + 5/18 * diag(ones(m - 1, 1), 1) + 4/9 * diag(ones(m, 1), 0); C3(1, 3) = 5/18; C3(m, m - 2) = 5/18; + C4 = 5/28 * diag(ones(m - 2, 1), -2) + 9/28 * diag(ones(m - 1, 1), -1) + 5/28 * diag(ones(m - 1, 1), 1) + 9/28 * diag(ones(m, 1), 0); C4(2, 4) = 5/28; C4(1, 3) = 5/28; C4(1, 4) = 9/28; C4(m, m - 2) = 5/28; C4(m, m - 3) = 5/28; + C5 = 1/7 * diag(ones(m - 2, 1), -2) + 1/7 * diag(ones(m - 2, 1), 2) + 8/35 * diag(ones(m - 1, 1), -1) + 8/35 * diag(ones(m - 1, 1), 1) + 9/35 * diag(ones(m, 1), 0); C5(2, 4) = 1/7; C5(1, 3) = 1/7; C5(1, 4) = 1/7; C5(1, 5) = 8/35; C5(2, 5) = 1/7; C5(m - 1, m - 4) = 1/7; C5(m, m - 4) = 8/35; C5(m, m - 3) = 1/7; + + C2 = sparse(diag(C2 * c)); + C3 = sparse(diag(C3 * c)); + C4 = sparse(diag(C4 * c)); + C5 = sparse(diag(C5 * c)); + + % Remainder term added to wide second derivative operator + R = (1/30735936 / h) * transpose(DD_12) * C1 * DD_12 + (1/6403320 / h) * transpose(DD_11) * C2 * DD_11 + (1/1293600 / h) * transpose(DD_10) * C3 * DD_10 + (1/249480 / h) * transpose(DD_9) * C4 * DD_9 + (1/44352 / h) * transpose(DD_8) * C5 * DD_8; + D2 = D1 * C1 * D1 - H \ R; + end + + % Wide stencil + function D2 = D2_wide(c) + % Here we add variable diffusion + C1 = sparse(diag(c)); + D2 = D1 * C1 * D1; + end + + switch options.stencil_width + case 'minimal' + D2 = @D2_fun_minimal; + case 'nonminimal' + D2 = @D2_fun_nonminimal; + case 'wide' + D2 = @D2_fun_wide; + otherwise + error('No option %s for stencil width', options.stencil_width) + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Artificial dissipation operator %%% + switch options.AD + case 'upwind' + % This is the choice that yield 11th order Upwind + DI = H \ (transpose(DD_6) * DD_6) * (-1/5544); + case 'op' + % This choice will preserve the order of the underlying + % Non-dissipative D1 SBP operator + DI = H \ (transpose(DD_7) * DD_7) * (-1 / (5 * 5544)); + % Notice that you can use any negative number instead of (-1/(5*5544)) + otherwise + error("Artificial dissipation options '%s' not implemented.", option.AD) + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%% +end \ No newline at end of file