Mercurial > repos > public > sbplib

function [H, HI, D1, D2, DI] = d2_noneq_variable_10(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 = 10;
    order = 10;

    %%%% Norm matrix %%%%%%%%
    P = zeros(BP, 1);
    P0 = 1.0000000000000e-01;
    P1 = 5.8980851260667e-01;
    P2 = 9.5666820955973e-01;
    P3 = 1.1500297411596e+00;
    P4 = 1.1232986993248e+00;
    P5 = 1.0123020150951e+00;
    P6 = 9.9877122702527e-01;
    P7 = 1.0000873322761e+00;
    P8 = 1.0000045540888e+00;
    P9 = 9.9999861455083e-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/1260, 5/504, -5/84, 5/21, -5/6, 0, 5/6, -5/21, 5/84, -5/504, 1/1260];
    d = repmat(d, N, 1);
    Q = spdiags(d, -order / 2:order / 2, N, N);

    % Boundaries
    Q0_0 = -5.0000000000000e-01;
    Q0_1 = 6.7548747038002e-01;
    Q0_2 = -2.6691978151546e-01;
    Q0_3 = 1.4438714982130e-01;
    Q0_4 = -7.7273673750760e-02;
    Q0_5 = 2.5570078343005e-02;
    Q0_6 = 4.2808774693299e-03;
    Q0_7 = -8.2902108933389e-03;
    Q0_8 = 3.2031176427908e-03;
    Q0_9 = -4.4502749689556e-04;
    Q0_10 = 0.0000000000000e+00;
    Q0_11 = 0.0000000000000e+00;
    Q0_12 = 0.0000000000000e+00;
    Q0_13 = 0.0000000000000e+00;
    Q0_14 = 0.0000000000000e+00;
    Q1_0 = -6.7548747038002e-01;
    Q1_1 = 0.0000000000000e+00;
    Q1_2 = 9.5146052715180e-01;
    Q1_3 = -4.2442349882626e-01;
    Q1_4 = 2.1538865145190e-01;
    Q1_5 = -7.1939778160350e-02;
    Q1_6 = -8.2539187832840e-03;
    Q1_7 = 1.9930661669090e-02;
    Q1_8 = -7.7433256989613e-03;
    Q1_9 = 1.0681515760869e-03;
    Q1_10 = 0.0000000000000e+00;
    Q1_11 = 0.0000000000000e+00;
    Q1_12 = 0.0000000000000e+00;
    Q1_13 = 0.0000000000000e+00;
    Q1_14 = 0.0000000000000e+00;
    Q2_0 = 2.6691978151546e-01;
    Q2_1 = -9.5146052715180e-01;
    Q2_2 = 0.0000000000000e+00;
    Q2_3 = 9.6073770842387e-01;
    Q2_4 = -3.9378595264609e-01;
    Q2_5 = 1.3302097358959e-01;
    Q2_6 = 8.1200458151489e-05;
    Q2_7 = -2.3849770528789e-02;
    Q2_8 = 9.6600442856829e-03;
    Q2_9 = -1.3234579460680e-03;
    Q2_10 = 0.0000000000000e+00;
    Q2_11 = 0.0000000000000e+00;
    Q2_12 = 0.0000000000000e+00;
    Q2_13 = 0.0000000000000e+00;
    Q2_14 = 0.0000000000000e+00;
    Q3_0 = -1.4438714982130e-01;
    Q3_1 = 4.2442349882626e-01;
    Q3_2 = -9.6073770842387e-01;
    Q3_3 = 0.0000000000000e+00;
    Q3_4 = 9.1551097634196e-01;
    Q3_5 = -2.8541713079648e-01;
    Q3_6 = 4.1398809121293e-02;
    Q3_7 = 1.7256059167927e-02;
    Q3_8 = -9.4349194803610e-03;
    Q3_9 = 1.3875650645663e-03;
    Q3_10 = 0.0000000000000e+00;
    Q3_11 = 0.0000000000000e+00;
    Q3_12 = 0.0000000000000e+00;
    Q3_13 = 0.0000000000000e+00;
    Q3_14 = 0.0000000000000e+00;
    Q4_0 = 7.7273673750760e-02;
    Q4_1 = -2.1538865145190e-01;
    Q4_2 = 3.9378595264609e-01;
    Q4_3 = -9.1551097634196e-01;
    Q4_4 = 0.0000000000000e+00;
    Q4_5 = 8.3519401865051e-01;
    Q4_6 = -2.0586492924974e-01;
    Q4_7 = 3.1230261235901e-02;
    Q4_8 = -2.0969453466651e-04;
    Q4_9 = -5.0965470499782e-04;
    Q4_10 = 0.0000000000000e+00;
    Q4_11 = 0.0000000000000e+00;
    Q4_12 = 0.0000000000000e+00;
    Q4_13 = 0.0000000000000e+00;
    Q4_14 = 0.0000000000000e+00;
    Q5_0 = -2.5570078343005e-02;
    Q5_1 = 7.1939778160350e-02;
    Q5_2 = -1.3302097358959e-01;
    Q5_3 = 2.8541713079648e-01;
    Q5_4 = -8.3519401865051e-01;
    Q5_5 = 0.0000000000000e+00;
    Q5_6 = 8.1046389580138e-01;
    Q5_7 = -2.1879194972141e-01;
    Q5_8 = 5.2977237804899e-02;
    Q5_9 = -9.0146730522360e-03;
    Q5_10 = 7.9365079365079e-04;
    Q5_11 = 0.0000000000000e+00;
    Q5_12 = 0.0000000000000e+00;
    Q5_13 = 0.0000000000000e+00;
    Q5_14 = 0.0000000000000e+00;
    Q6_0 = -4.2808774693299e-03;
    Q6_1 = 8.2539187832840e-03;
    Q6_2 = -8.1200458151489e-05;
    Q6_3 = -4.1398809121293e-02;
    Q6_4 = 2.0586492924974e-01;
    Q6_5 = -8.1046389580138e-01;
    Q6_6 = 0.0000000000000e+00;
    Q6_7 = 8.2787884456005e-01;
    Q6_8 = -2.3582460382545e-01;
    Q6_9 = 5.9178678209520e-02;
    Q6_10 = -9.9206349206349e-03;
    Q6_11 = 7.9365079365079e-04;
    Q6_12 = 0.0000000000000e+00;
    Q6_13 = 0.0000000000000e+00;
    Q6_14 = 0.0000000000000e+00;
    Q7_0 = 8.2902108933389e-03;
    Q7_1 = -1.9930661669090e-02;
    Q7_2 = 2.3849770528789e-02;
    Q7_3 = -1.7256059167927e-02;
    Q7_4 = -3.1230261235901e-02;
    Q7_5 = 2.1879194972141e-01;
    Q7_6 = -8.2787884456005e-01;
    Q7_7 = 0.0000000000000e+00;
    Q7_8 = 8.3301028859275e-01;
    Q7_9 = -2.3804321850015e-01;
    Q7_10 = 5.9523809523809e-02;
    Q7_11 = -9.9206349206349e-03;
    Q7_12 = 7.9365079365079e-04;
    Q7_13 = 0.0000000000000e+00;
    Q7_14 = 0.0000000000000e+00;
    Q8_0 = -3.2031176427908e-03;
    Q8_1 = 7.7433256989613e-03;
    Q8_2 = -9.6600442856829e-03;
    Q8_3 = 9.4349194803610e-03;
    Q8_4 = 2.0969453466651e-04;
    Q8_5 = -5.2977237804899e-02;
    Q8_6 = 2.3582460382545e-01;
    Q8_7 = -8.3301028859275e-01;
    Q8_8 = 0.0000000000000e+00;
    Q8_9 = 8.3333655748509e-01;
    Q8_10 = -2.3809523809524e-01;
    Q8_11 = 5.9523809523809e-02;
    Q8_12 = -9.9206349206349e-03;
    Q8_13 = 7.9365079365079e-04;
    Q8_14 = 0.0000000000000e+00;
    Q9_0 = 4.4502749689556e-04;
    Q9_1 = -1.0681515760869e-03;
    Q9_2 = 1.3234579460680e-03;
    Q9_3 = -1.3875650645663e-03;
    Q9_4 = 5.0965470499782e-04;
    Q9_5 = 9.0146730522360e-03;
    Q9_6 = -5.9178678209520e-02;
    Q9_7 = 2.3804321850015e-01;
    Q9_8 = -8.3333655748509e-01;
    Q9_9 = 0.0000000000000e+00;
    Q9_10 = 8.3333333333333e-01;
    Q9_11 = -2.3809523809524e-01;
    Q9_12 = 5.9523809523809e-02;
    Q9_13 = -9.9206349206349e-03;
    Q9_14 = 7.9365079365079e-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_5 = (-diag(ones(m - 3, 1), -3) + 5 * diag(ones(m - 2, 1), -2) - 10 * diag(ones(m - 1, 1), -1) + 10 * diag(ones(m, 1), 0) - 5 * diag(ones(m - 1, 1), 1) + diag(ones(m - 2, 1), 2));
    DD_5(1:8, 1:10) = [0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; -0.88425676349344827581e1 0.18735837038089672208e2 -0.17076524806228533350e2 0.10664693462778093925e2 -0.43403898038449865885e1 0.85895174414023656383e0 0 0 0 0; 0 -0.13262411219221127021e1 0.45552739083119183335e1 -0.73424543813164879330e1 0.70876331528332015948e1 -0.38059884697709954324e1 0.83177691186447613913e0 0 0 0; 0 0 -0.67600766166626540709e0 0.32310531977344450754e1 -0.69639522132755183916e1 0.77488870404680024829e1 -0.42187263911386203939e1 0.87874602787795663432e0 0 0; 0 0 0 -0.66326118352702166971e0 0.38132489024062224761e1 -0.84909798376351149047e1 0.90434791287189283383e1 -0.46461964978830115989e1 0.94370948791999735890e0 0; 0 0 0 0 -0.86903681018048694197e0 0.47050202961852482685e1 -0.96960545214617434798e1 0.97952957162584740145e1 -0.49228410242754295559e1 0.98761634347393769455e0; ];
    DD_5(m - 6:m, m - 9:m) = [-0.98761634347393769455e0 0.49228410242754295559e1 -0.97952957162584740145e1 0.96960545214617434798e1 -0.47050202961852482685e1 0.86903681018048694197e0 0 0 0 0; 0 -0.94370948791999735890e0 0.46461964978830115989e1 -0.90434791287189283383e1 0.84909798376351149047e1 -0.38132489024062224761e1 0.66326118352702166971e0 0 0 0; 0 0 -0.87874602787795663432e0 0.42187263911386203939e1 -0.77488870404680024829e1 0.69639522132755183916e1 -0.32310531977344450754e1 0.67600766166626540709e0 0 0; 0 0 0 -0.83177691186447613913e0 0.38059884697709954324e1 -0.70876331528332015948e1 0.73424543813164879330e1 -0.45552739083119183335e1 0.13262411219221127021e1 0; 0 0 0 0 -0.85895174414023656383e0 0.43403898038449865885e1 -0.10664693462778093925e2 0.17076524806228533350e2 -0.18735837038089672208e2 0.88425676349344827581e1; 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0; ];
    DD_5 = sparse(DD_5);

    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:8, 1:11) = [0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; 0.97690498198305661286e1 -0.22164087301995739032e2 0.23898275711233304083e2 -0.19893851160044639280e2 0.12625396854743289626e2 -0.51537104648414193830e1 0.91892654107463785790e0 0 0 0 0; 0 0.13105269062480722931e1 -0.51692977530964314903e1 0.10448225399376470222e2 -0.13885092532071720368e2 0.11417965409312986297e2 -0.49906614711868568348e1 0.86833404141747988109e0 0 0 0; 0 0 0.64511698871199831925e0 -0.37163607887886698482e1 0.10284728893981919286e2 -0.15497774080936004966e2 0.12656179173415861182e2 -0.52724761672677398059e1 0.90058598088263583338e0 0 0; 0 0 0 0.64016413450696574250e0 -0.45192323253005276829e1 0.12736469756452672357e2 -0.18086958257437856677e2 0.13938589493649034797e2 -0.56622569275199841534e1 0.95322412564969561677e0 0; 0 0 0 0 0.86005005088560527649e0 -0.56460243554222979222e1 0.14544081782192615220e2 -0.19590591432516948029e2 0.14768523072826288668e2 -0.59256980608436261673e1 0.98965894287836295484e0; ];
    DD_6(m - 7:m, m - 10:m) = [0.98965894287836295484e0 -0.59256980608436261673e1 0.14768523072826288668e2 -0.19590591432516948029e2 0.14544081782192615220e2 -0.56460243554222979222e1 0.86005005088560527649e0 0 0 0 0; 0 0.95322412564969561677e0 -0.56622569275199841534e1 0.13938589493649034797e2 -0.18086958257437856677e2 0.12736469756452672357e2 -0.45192323253005276829e1 0.64016413450696574250e0 0 0 0; 0 0 0.90058598088263583338e0 -0.52724761672677398059e1 0.12656179173415861182e2 -0.15497774080936004966e2 0.10284728893981919286e2 -0.37163607887886698482e1 0.64511698871199831925e0 0 0; 0 0 0 0.86833404141747988109e0 -0.49906614711868568348e1 0.11417965409312986297e2 -0.13885092532071720368e2 0.10448225399376470222e2 -0.51692977530964314903e1 0.13105269062480722931e1 0; 0 0 0 0 0.91892654107463785790e0 -0.51537104648414193830e1 0.12625396854743289626e2 -0.19893851160044639280e2 0.23898275711233304083e2 -0.22164087301995739032e2 0.97690498198305661286e1; 0 0 0 0 0 0 0 0 0 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: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 0 0 0 0 0 0 0 0 0 0 0; -0.10633444157744515342e2 0.25551717302255079868e2 -0.31639558087487298930e2 0.33026832975062979647e2 -0.28856215669710369771e2 0.18037986626944967840e2 -0.64324857875224650053e1 0.94516679820162169311e0 0 0 0 0; 0 -0.12971946051930944772e1 0.57552633215463132945e1 -0.14020486493241325087e2 0.23923936099960072962e2 -0.26641919288396968027e2 0.17467315149153998922e2 -0.60783382899223591677e1 0.89142410609336157976e0 0 0 0; 0 0 -0.61968315701047317846e0 0.41847682905586435755e1 -0.14220312881453249947e2 0.27121104641638008690e2 -0.29531084737970342758e2 0.18453666585437089321e2 -0.63041018661784508337e1 0.91564312497877513070e0 0 0; 0 0 0 -0.62096054671195295664e0 0.52179151332917540372e1 -0.17831057659033741300e2 0.31652176950516249184e2 -0.32523375485181081192e2 0.19817899246319944537e2 -0.66725688795478693174e1 0.95997124034669700791e0 0; 0 0 0 0 -0.85241549236735026150e0 0.65870284146593475759e1 -0.20361714495069661308e2 0.34283535006904659051e2 -0.34459887169928006891e2 0.20739943212952691585e2 -0.69276126001485406839e1 0.99112312299686093167e0; ];
    DD_7(m - 7:m, m - 11:m) = [-0.99112312299686093167e0 0.69276126001485406839e1 -0.20739943212952691585e2 0.34459887169928006891e2 -0.34283535006904659051e2 0.20361714495069661308e2 -0.65870284146593475759e1 0.85241549236735026150e0 0 0 0 0; 0 -0.95997124034669700791e0 0.66725688795478693174e1 -0.19817899246319944537e2 0.32523375485181081192e2 -0.31652176950516249184e2 0.17831057659033741300e2 -0.52179151332917540372e1 0.62096054671195295664e0 0 0 0; 0 0 -0.91564312497877513070e0 0.63041018661784508337e1 -0.18453666585437089321e2 0.29531084737970342758e2 -0.27121104641638008690e2 0.14220312881453249947e2 -0.41847682905586435755e1 0.61968315701047317846e0 0 0; 0 0 0 -0.89142410609336157976e0 0.60783382899223591677e1 -0.17467315149153998922e2 0.26641919288396968027e2 -0.23923936099960072962e2 0.14020486493241325087e2 -0.57552633215463132945e1 0.12971946051930944772e1 0; 0 0 0 0 -0.94516679820162169311e0 0.64324857875224650053e1 -0.18037986626944967840e2 0.28856215669710369771e2 -0.33026832975062979647e2 0.31639558087487298930e2 -0.25551717302255079868e2 0.10633444157744515342e2; 0 0 0 0 0 0 0 0 0 0 0 0; 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:9, 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.11447706798618566234e2 -0.28904884139025747734e2 0.40258353260409964223e2 -0.50649997399180126990e2 0.56821824921675878869e2 -0.48101297671853247575e2 0.25729943150089860021e2 -0.75613343856129735449e1 0.95968546487782649665e0 0 0 0 0; 0 0.12856328177474934636e1 -0.63181271116961104382e1 0.18042992864608611563e2 -0.37804272468074727719e2 0.53283838576793936053e2 -0.46579507064410663791e2 0.24313353159689436671e2 -0.71313928487468926381e1 0.90748207408891683663e0 0 0 0; 0 0 0.59820007352804876260e0 -0.46391245450012267365e1 0.18764346418211454409e2 -0.43393767426620813904e2 0.59062169475940685515e2 -0.49209777561165571522e2 0.25216407464713803335e2 -0.73251449998302010456e1 0.92669110022382118702e0 0 0; 0 0 0 0.60460011916248255841e0 -0.59103958199322344450e1 0.23774743545378321733e2 -0.50643483120825998695e2 0.65046750970362162385e2 -0.52847731323519852099e2 0.26690275518191477270e2 -0.76797699227735760633e1 0.96501003395721735560e0 0; 0 0 0 0 0.84578719850252712660e0 -0.75280324738963972296e1 0.27148952660092881743e2 -0.54853656011047454481e2 0.68919774339856013782e2 -0.55306515234540510895e2 0.27710450400594162736e2 -0.79289849839748874534e1 0.99222410441366467138e0; ];
    DD_8(m - 8:m, m - 12:m) = [0.99222410441366467138e0 -0.79289849839748874534e1 0.27710450400594162736e2 -0.55306515234540510895e2 0.68919774339856013782e2 -0.54853656011047454481e2 0.27148952660092881743e2 -0.75280324738963972296e1 0.84578719850252712660e0 0 0 0 0; 0 0.96501003395721735560e0 -0.76797699227735760633e1 0.26690275518191477270e2 -0.52847731323519852099e2 0.65046750970362162385e2 -0.50643483120825998695e2 0.23774743545378321733e2 -0.59103958199322344450e1 0.60460011916248255841e0 0 0 0; 0 0 0.92669110022382118702e0 -0.73251449998302010456e1 0.25216407464713803335e2 -0.49209777561165571522e2 0.59062169475940685515e2 -0.43393767426620813904e2 0.18764346418211454409e2 -0.46391245450012267365e1 0.59820007352804876260e0 0 0; 0 0 0 0.90748207408891683663e0 -0.71313928487468926381e1 0.24313353159689436671e2 -0.46579507064410663791e2 0.53283838576793936053e2 -0.37804272468074727719e2 0.18042992864608611563e2 -0.63181271116961104382e1 0.12856328177474934636e1 0; 0 0 0 0 0.95968546487782649665e0 -0.75613343856129735449e1 0.25729943150089860021e2 -0.48101297671853247575e2 0.56821824921675878869e2 -0.50649997399180126990e2 0.40258353260409964223e2 -0.28904884139025747734e2 0.11447706798618566234e2; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 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: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 0 0 0 0 0 0 0 0 0 0 0 0 0; -0.12220346479758703402e2 0.32228164479135716142e2 -0.49720065159556828467e2 0.73329283892516285142e2 -0.10101269332559123572e3 0.10822791976166980704e3 -0.77189829450269580063e2 0.34026004735258380952e2 -0.86371691839004384699e1 0.96873073049659684460e0 0 0 0 0; 0 -0.12754371756934054072e1 0.68614776237194576426e1 -0.22502238094968228382e2 0.56120004490560654424e2 -0.95910909438229084895e2 0.10480389089492399353e3 -0.72940059479068310012e2 0.32091267819361016871e2 -0.81673386668002515297e1 0.91934202619415775756e0 0 0 0; 0 0 -0.57969473693003825890e0 0.50815098898235167347e1 -0.23911428372444897689e2 0.65090651139931220857e2 -0.10631190505669323393e3 0.11072199951262253592e3 -0.75649222394141410004e2 0.32963152499235904705e2 -0.83402199020143906832e1 0.93515742061079234147e0 0 0; 0 0 0 -0.59039909771982400264e0 0.65974918490578446834e1 -0.30567527415486413657e2 0.75965224681238998042e2 -0.11708415174665189229e3 0.11890739547791966722e3 -0.80070826554574431809e2 0.34558964652481092285e2 -0.86850903056149562004e1 0.96891845934991572940e0 0; 0 0 0 0 -0.83993614063969059233e0 0.84690365331334468834e1 -0.34905796277262276527e2 0.82280484016571181722e2 -0.12405559381174082481e3 0.12443965927771614951e3 -0.83131351201782488207e2 0.35680432427886993540e2 -0.89300169397229820424e1 0.99308211584049051802e0; ];
    DD_9(m - 8:m, m - 13:m) = [-0.99308211584049051802e0 0.89300169397229820424e1 -0.35680432427886993540e2 0.83131351201782488207e2 -0.12443965927771614951e3 0.12405559381174082481e3 -0.82280484016571181722e2 0.34905796277262276527e2 -0.84690365331334468834e1 0.83993614063969059233e0 0 0 0 0; 0 -0.96891845934991572940e0 0.86850903056149562004e1 -0.34558964652481092285e2 0.80070826554574431809e2 -0.11890739547791966722e3 0.11708415174665189229e3 -0.75965224681238998042e2 0.30567527415486413657e2 -0.65974918490578446834e1 0.59039909771982400264e0 0 0 0; 0 0 -0.93515742061079234147e0 0.83402199020143906832e1 -0.32963152499235904705e2 0.75649222394141410004e2 -0.11072199951262253592e3 0.10631190505669323393e3 -0.65090651139931220857e2 0.23911428372444897689e2 -0.50815098898235167347e1 0.57969473693003825890e0 0 0; 0 0 0 -0.91934202619415775756e0 0.81673386668002515297e1 -0.32091267819361016871e2 0.72940059479068310012e2 -0.10480389089492399353e3 0.95910909438229084895e2 -0.56120004490560654424e2 0.22502238094968228382e2 -0.68614776237194576426e1 0.12754371756934054072e1 0; 0 0 0 0 -0.96873073049659684460e0 0.86371691839004384699e1 -0.34026004735258380952e2 0.77189829450269580063e2 -0.10822791976166980704e3 0.10101269332559123572e3 -0.73329283892516285142e2 0.49720065159556828467e2 -0.32228164479135716142e2 0.12220346479758703402e2; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 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:10, 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.12957678688124699986e2 -0.35525089722887367966e2 0.59995471673283665865e2 -0.10161365491269611714e3 0.16661352548983481943e3 -0.21645583952333961409e3 0.19297457362567395016e3 -0.11342001578419460317e3 0.43185845919502192349e2 -0.96873073049659684460e1 0.97481185166434302666e0 0 0 0 0; 0 0.12663266431786644877e1 -0.73880195719183439765e1 0.27386715439651102593e2 -0.79459763018974764721e2 0.15985151573038180816e3 -0.20960778178984798706e3 0.18235014869767077503e3 -0.10697089273120338957e3 0.40836693334001257648e2 -0.91934202619415775756e1 0.92847752900245498778e0 0 0 0; 0 0 0.56350506801536116666e0 -0.55135043604652958562e1 0.29656869502625787574e2 -0.92986644485616029795e2 0.17718650842782205655e3 -0.22144399902524507185e3 0.18912305598535352501e3 -0.10987717499745301568e3 0.41701099510071953416e2 -0.93515742061079234147e1 0.94185858099865288757e0 0 0; 0 0 0 0.57788899624281661128e0 -0.72798346274475049971e1 0.38209409269358017071e2 -0.10852174954462714006e3 0.19514025291108648715e3 -0.23781479095583933444e3 0.20017706638643607952e3 -0.11519654884160364095e3 0.43425451528074781002e2 -0.96891845934991572940e1 0.97203947181859638306e0 0; 0 0 0 0 0.83470299758184640199e0 -0.94100405923704965371e1 0.43632245346577845659e2 -0.11754354859510168817e3 0.20675932301956804135e3 -0.24887931855543229903e3 0.20782837800445622052e3 -0.11893477475962331180e3 0.44650084698614910212e2 -0.99308211584049051802e1 0.99376959413383658153e0; ];
    DD_10(m - 9:m, m - 14:m) = [0.99376959413383658153e0 -0.99308211584049051802e1 0.44650084698614910212e2 -0.11893477475962331180e3 0.20782837800445622052e3 -0.24887931855543229903e3 0.20675932301956804135e3 -0.11754354859510168817e3 0.43632245346577845659e2 -0.94100405923704965371e1 0.83470299758184640199e0 0 0 0 0; 0 0.97203947181859638306e0 -0.96891845934991572940e1 0.43425451528074781002e2 -0.11519654884160364095e3 0.20017706638643607952e3 -0.23781479095583933444e3 0.19514025291108648715e3 -0.10852174954462714006e3 0.38209409269358017071e2 -0.72798346274475049971e1 0.57788899624281661128e0 0 0 0; 0 0 0.94185858099865288757e0 -0.93515742061079234147e1 0.41701099510071953416e2 -0.10987717499745301568e3 0.18912305598535352501e3 -0.22144399902524507185e3 0.17718650842782205655e3 -0.92986644485616029795e2 0.29656869502625787574e2 -0.55135043604652958562e1 0.56350506801536116666e0 0 0; 0 0 0 0.92847752900245498778e0 -0.91934202619415775756e1 0.40836693334001257648e2 -0.10697089273120338957e3 0.18235014869767077503e3 -0.20960778178984798706e3 0.15985151573038180816e3 -0.79459763018974764721e2 0.27386715439651102593e2 -0.73880195719183439765e1 0.12663266431786644877e1 0; 0 0 0 0 0.97481185166434302666e0 -0.96873073049659684460e1 0.43185845919502192349e2 -0.11342001578419460317e3 0.19297457362567395016e3 -0.21645583952333961409e3 0.16661352548983481943e3 -0.10161365491269611714e3 0.59995471673283665865e2 -0.35525089722887367966e2 0.12957678688124699986e2; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 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);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%% Difference operator %%
    D1 = H \ Q;

    % Helper functions for constructing D2(c)
    % TODO: Consider changing sparse(diag(...)) to spdiags(....)

    % Minimal 11 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 = 2/7 * diag(ones(m - 1, 1), -1) + 2/7 * diag(ones(m - 1, 1), 1) + 3/7 * diag(ones(m, 1), 0); C3(1, 3) = 2/7; C3(m, m - 2) = 2/7;
        C4 = 1/5 * diag(ones(m - 2, 1), -2) + 3/10 * diag(ones(m - 1, 1), -1) + 1/5 * diag(ones(m - 1, 1), 1) + 3/10 * diag(ones(m, 1), 0); C4(2, 4) = 1/5; C4(1, 3) = 1/5; C4(1, 4) = 3/10; C4(m, m - 2) = 1/5; C4(m, m - 3) = 1/5;
        C5 = 1/5 * diag(ones(m - 2, 1), -2) + 1/5 * diag(ones(m - 2, 1), 2) + 1/5 * diag(ones(m - 1, 1), -1) + 1/5 * diag(ones(m - 1, 1), 1) + 1/5 * diag(ones(m, 1), 0); C5(2, 4) = 1/5; C5(1, 4) = 1/5; C5(1, 5) = 1/5; C5(2, 5) = 1/5; C5(m - 1, m - 4) = 1/5; C5(m, m - 4) = 1/5; C5(m, m - 3) = 1/5;

        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/1587600 / h) * transpose(DD_10) * C1 * DD_10 + (1/317520 / h) * transpose(DD_9) * C2 * DD_9 + (1/60480 / h) * transpose(DD_8) * C3 * DD_8 + (1/10584 / h) * transpose(DD_7) * C4 * DD_7 + (1/1512 / h) * transpose(DD_6) * C5 * DD_6;
        D2 = D1 * C1 * D1 - H \ R;
    end

    % Few additional grid point in interior stencil cmp. to minimal
    function D2 = D2_fun_nonminimal(c)
        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 = 2/7 * diag(ones(m - 1, 1), -1) + 2/7 * diag(ones(m - 1, 1), 1) + 3/7 * diag(ones(m, 1), 0); C3(1, 3) = 2/7; C3(m, m - 2) = 2/7;
        C4 = 1/5 * diag(ones(m - 2, 1), -2) + 3/10 * diag(ones(m - 1, 1), -1) + 1/5 * diag(ones(m - 1, 1), 1) + 3/10 * diag(ones(m, 1), 0); C4(2, 4) = 1/5; C4(1, 3) = 1/5; C4(1, 4) = 3/10; C4(m, m - 2) = 1/5; C4(m, m - 3) = 1/5;

        C2 = sparse(diag(C2 * c));
        C3 = sparse(diag(C3 * c));
        C4 = sparse(diag(C4 * c));

        % Remainder term added to wide second derivative operator
        R = (1/1587600 / h) * transpose(DD_10) * C1 * DD_10 + (1/317520 / h) * transpose(DD_9) * C2 * DD_9 + (1/60480 / h) * transpose(DD_8) * C3 * DD_8 + (1/10584 / h) * transpose(DD_7) * C4 * DD_7;
        D2 = D1 * C1 * D1 - H \ R;
    end

    % Wide stencil
    function D2 = D2_fun_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 9th order Upwind
            DI = H \ (transpose(DD_5) * DD_5) * (-1/1260);
        case 'op'
            % This choice will preserve the order of the underlying
            % Non-dissipative D1 SBP operator
            DI = H \ (transpose(DD_6) * DD_6) * (-1 / (5 * 1260));
            % Notice that you can use any negative number instead of (-1/(5*1260))
        otherwise
            error("Artificial dissipation options '%s' not implemented.", option.AD)
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%
end
author	Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date	Sun, 13 Feb 2022 20:58:35 +0100
parents	1b0f2415237f
children	855871e0b852