Mercurial > repos > public > sbplib

view +sbp/+implementations/d2_blocknorm_8.m @ 267:f7ac3cd6eeaa operator_remake

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Sparsified all implementation files, removed all matlab warnings, fixed small bugs on minimum grid points.
author Martin Almquist <martin.almquist@it.uu.se>
date Fri, 09 Sep 2016 14:53:41 +0200
parents bfa130b7abf6
children
line wrap: on
line source

function [H, HI, D1, D2, e_1, e_m, M, Q, S_1, S_m] = d2_blocknorm_8(m,h)
    % Eighth order
    
    BP = 8;
    if(m<2*BP)
        error(['Operator requires at least ' num2str(2*BP) ' grid points']);
    end
    
    e_1=sparse(m,1);e_1(1)=1;
    e_m=sparse(m,1);e_m(m)=1;

    H_U=[     0.704266523e9 / 0.4180377600e10     0.4579586639e10 / 0.16257024000e11     -0.3623870581e10 / 0.9754214400e10    0.12753127559e11 / 0.29262643200e11    -0.3687413731e10 / 0.9754214400e10         0.2169892891e10 / 0.9754214400e10     -0.13224544841e11 / 0.146313216000e12        0.4142047e7 / 0.199065600e9;
            0.4579586639e10 / 0.16257024000e11   0.36543258551e11 / 0.20901888000e11     -0.8235820121e10 / 0.6967296000e10     0.1800520829e10 / 0.1393459200e10     -0.3725834681e10 / 0.4180377600e10         0.2588501879e10 / 0.6967296000e10         -0.10477621e8 / 0.995328000e9       -0.6589395529e10 / 0.146313216000e12;
           -0.3623870581e10 / 0.9754214400e10    -0.8235820121e10 / 0.6967296000e10         0.765685439e9 / 0.258048000e9      -0.3254203513e10 / 0.1393459200e10      0.2477267447e10 / 0.1393459200e10             -0.456533e6 / 0.491520e6              0.873043831e9 / 0.6967296000e10      0.4279558279e10 / 0.48771072000e11;
           0.12753127559e11 / 0.29262643200e11    0.1800520829e10 / 0.1393459200e10      -0.3254203513e10 / 0.1393459200e10    0.16428690611e11 / 0.4180377600e10       -0.460605929e9 / 0.199065600e9           0.1897042423e10 / 0.1393459200e10      -0.1151273401e10 / 0.4180377600e10        -0.65906413e8 / 0.650280960e9;
           -0.3687413731e10 / 0.9754214400e10    -0.3725834681e10 / 0.4180377600e10       0.2477267447e10 / 0.1393459200e10      -0.460605929e9 / 0.199065600e9       0.11663916211e11 / 0.4180377600e10          -0.307273957e9 / 0.278691840e9           0.117995903e9 / 0.464486400e9           0.95035807e8 / 0.1170505728e10;
            0.2169892891e10 / 0.9754214400e10     0.2588501879e10 / 0.6967296000e10           -0.456533e6 / 0.491520e6          0.1897042423e10 / 0.1393459200e10       -0.307273957e9 / 0.278691840e9              0.17436823e8 / 0.10321920e8         -0.1274455129e10 / 0.6967296000e10       -0.338917493e9 / 0.9754214400e10;
          -0.13224544841e11 / 0.146313216000e12     -0.10477621e8 / 0.995328000e9           0.873043831e9 / 0.6967296000e10    -0.1151273401e10 / 0.4180377600e10        0.117995903e9 / 0.464486400e9          -0.1274455129e10 / 0.6967296000e10      0.22041718711e11 / 0.20901888000e11       0.468461293e9 / 0.48771072000e11;
                0.4142047e7 / 0.199065600e9      -0.6589395529e10 / 0.146313216000e12     0.4279558279e10 / 0.48771072000e11      -0.65906413e8 / 0.650280960e9           0.95035807e8 / 0.1170505728e10          -0.338917493e9 / 0.9754214400e10         0.468461293e9 / 0.48771072000e11    0.20832744839e11 / 0.20901888000e11;
        ];


    H=speye(m);
    H(1:8,1:8)=H_U;
    H(m-7:m,m-7:m)=rot90( H_U(1:8,1:8) ,2 );
    H=H*h;
    HI=inv(H);

%     Q=-(1/280*diag(ones(m-4,1),4)-4/105*diag(ones(m-3,1),3)+1/5*diag(ones(m-2,1),2)-4/5*diag(ones(m-1,1),1)+4/5*diag(ones(m-1,1),-1)-1/5*diag(ones(m-2,1),-2)+4/105*diag(ones(m-3,1),-3)-1/280*diag(ones(m-4,1),-4));
    diags   = -4:4;
    stencil = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280];
    Q = stripeMatrix(stencil, diags, m);

    Q_U = [-0.1e1 / 0.2e1 0.16262381e8 / 0.15876000e8 -0.3770744693e10 / 0.3048192000e10 0.290431859e9 / 0.203212800e9 -0.363704879e9 / 0.304819200e9 0.206906927e9 / 0.304819200e9 -0.248883679e9 / 0.1016064000e10 0.2665637e7 / 0.62208000e8; -0.16262381e8 / 0.15876000e8 0 0.138589901e9 / 0.62208000e8 -0.31764881e8 / 0.12441600e8 0.461249e6 / 0.193536e6 -0.347232997e9 / 0.217728000e9 0.1788157e7 / 0.2488320e7 -0.5926349e7 / 0.37632000e8; 0.3770744693e10 / 0.3048192000e10 -0.138589901e9 / 0.62208000e8 0 0.11741773e8 / 0.5443200e7 -0.39109817e8 / 0.17418240e8 0.10216441e8 / 0.5376000e7 -0.245131109e9 / 0.217728000e9 0.92809903e8 / 0.304819200e9; -0.290431859e9 / 0.203212800e9 0.31764881e8 / 0.12441600e8 -0.11741773e8 / 0.5443200e7 0 0.4634999e7 / 0.2488320e7 -0.144219869e9 / 0.87091200e8 0.17445643e8 / 0.14515200e8 -0.38142949e8 / 0.101606400e9; 0.363704879e9 / 0.304819200e9 -0.461249e6 / 0.193536e6 0.39109817e8 / 0.17418240e8 -0.4634999e7 / 0.2488320e7 0 0.7992221e7 / 0.5443200e7 -0.817951e6 / 0.829440e6 0.4455517e7 / 0.13547520e8; -0.206906927e9 / 0.304819200e9 0.347232997e9 / 0.217728000e9 -0.10216441e8 / 0.5376000e7 0.144219869e9 / 0.87091200e8 -0.7992221e7 / 0.5443200e7 0 0.68487373e8 / 0.62208000e8 -0.1032638773e10 / 0.3048192000e10; 0.248883679e9 / 0.1016064000e10 -0.1788157e7 / 0.2488320e7 0.245131109e9 / 0.217728000e9 -0.17445643e8 / 0.14515200e8 0.817951e6 / 0.829440e6 -0.68487373e8 / 0.62208000e8 0 0.39529771e8 / 0.47628000e8; -0.2665637e7 / 0.62208000e8 0.5926349e7 / 0.37632000e8 -0.92809903e8 / 0.304819200e9 0.38142949e8 / 0.101606400e9 -0.4455517e7 / 0.13547520e8 0.1032638773e10 / 0.3048192000e10 -0.39529771e8 / 0.47628000e8 0;];

    Q(1:8,1:8)=Q_U;
    Q(m-7:m,m-7:m)=rot90( -Q_U(1:8,1:8) ,2 );

    D1=H\Q;


    M_U =[0.27667249117e11 / 0.18289152000e11 -0.17100791927e11 / 0.6096384000e10 0.6123596021e10 / 0.2032128000e10 -0.12420079921e11 / 0.3657830400e10 0.3352522937e10 / 0.1219276800e10 -0.3030351383e10 / 0.2032128000e10 0.8955233071e10 / 0.18289152000e11 -0.448917533e9 / 0.6096384000e10; -0.2443029521e10 / 0.870912000e9 0.3279926909e10 / 0.373248000e9 -0.150833107e9 / 0.11612160e8 0.2418903029e10 / 0.174182400e9 -0.1195687489e10 / 0.104509440e9 0.1864443097e10 / 0.290304000e9 -0.275412413e9 / 0.124416000e9 0.26267539e8 / 0.74649600e8; 0.875033123e9 / 0.290304000e9 -0.754432799e9 / 0.58060800e8 0.262316881e9 / 0.10752000e8 -0.1615952663e10 / 0.58060800e8 0.1304948581e10 / 0.58060800e8 -0.59605951e8 / 0.4608000e7 0.270029509e9 / 0.58060800e8 -0.225377137e9 / 0.290304000e9; -0.71086111e8 / 0.20901888e8 0.2425994741e10 / 0.174182400e9 -0.1622486807e10 / 0.58060800e8 0.735382895e9 / 0.20901888e8 -0.1016121419e10 / 0.34836480e8 0.190014817e9 / 0.11612160e8 -0.447155539e9 / 0.74649600e8 0.179406911e9 / 0.174182400e9; 0.480578879e9 / 0.174182400e9 -0.6018333509e10 / 0.522547200e9 0.1319413093e10 / 0.58060800e8 -0.205032463e9 / 0.6967296e7 0.551889007e9 / 0.20901888e8 -0.887809303e9 / 0.58060800e8 0.914606453e9 / 0.174182400e9 -0.67482881e8 / 0.74649600e8; -0.434493809e9 / 0.290304000e9 0.1878773977e10 / 0.290304000e9 -0.423185977e9 / 0.32256000e8 0.964538597e9 / 0.58060800e8 -0.894343447e9 / 0.58060800e8 0.345461491e9 / 0.32256000e8 -0.1288081307e10 / 0.290304000e9 0.199200163e9 / 0.290304000e9; 0.183060319e9 / 0.373248000e9 -0.276656573e9 / 0.124416000e9 0.54579137e8 / 0.11612160e8 -0.3169984837e10 / 0.522547200e9 0.184339633e9 / 0.34836480e8 -0.1289417627e10 / 0.290304000e9 0.1431981949e10 / 0.373248000e9 -0.307164061e9 / 0.174182400e9; -0.449332253e9 / 0.6096384000e10 0.1290053923e10 / 0.3657830400e10 -0.1588745239e10 / 0.2032128000e10 0.1267377593e10 / 0.1219276800e10 -0.3326650673e10 / 0.3657830400e10 0.1396036981e10 / 0.2032128000e10 -0.2150231371e10 / 0.1219276800e10 0.52544801501e11 / 0.18289152000e11;];

%     M=-(-1/560*diag(ones(m-4,1),4)+8/315*diag(ones(m-3,1),3)-1/5*diag(ones(m-2,1),2)+8/5*diag(ones(m-1,1),1)+8/5*diag(ones(m-1,1),-1)-1/5*diag(ones(m-2,1),-2)+8/315*diag(ones(m-3,1),-3)-1/560*diag(ones(m-4,1),-4)-205/72*diag(ones(m,1),0));
    diags   = -4:4;
    left_stencil = [-1/560,8/315,-1/5,8/5];
    stencil = [left_stencil,-205/72,fliplr(left_stencil)];
    M = stripeMatrix(stencil, diags, m);

    M(1:8,1:8)=M_U;

    M(m-7:m,m-7:m)=rot90(  M_U ,2 );
    M=M/h;

    % DS_U=[0.363e3 / 0.140e3 -7 0.21e2 / 0.2e1 -0.35e2 / 0.3e1 0.35e2 / 0.4e1 -0.21e2 / 0.5e1 0.7e1 / 0.6e1 -0.1e1 / 0.7e1;];
    % DS=sparse(m,m);
    % DS(1,1:8)=DS_U;
    % DS(m,m-7:m)=fliplr(DS_U);
    % DS=DS/h;

    %D2=HI*(-M+DS);

    S_U=-[0.363e3 / 0.140e3 -7 0.21e2 / 0.2e1 -0.35e2 / 0.3e1 0.35e2 / 0.4e1 -0.21e2 / 0.5e1 0.7e1 / 0.6e1 -0.1e1 / 0.7e1;]/h;
    S_1=sparse(1,m);
    S_1(1:8)=S_U;
    S_m=sparse(1,m);
    S_m(m-7:m)=fliplr(-S_U);


    D2=H\(-M - e_1*S_1+e_m*S_m);
    S_1 = S_1';
    S_m = S_m';
end