view +sbp/ordinary12.m @ 255:df3cc9c5dffc operator_remake

Added ordinary 12th order accurate, with D1*D1 as 2nd derivative.
author Martin Almquist <martin.almquist@it.uu.se>
date Wed, 07 Sep 2016 15:54:41 +0200
parents
children 5714fda0db85
line wrap: on
line source

function [H, HI, D1, D2, e_1, e_m, M, Q, S_1, S_m] = ordinary12(m,h)

% D2 = D1*D1, wide!

H=diag(ones(m,1),0);
H(1:15,1:15)=diag([2.880607858916397e-01,...
    1.560376162339675e+00,...
    2.403139445479289e-01,...
    1.531292896754208e+00,...
    1.245107365551116e+00,...
    2.403139445479289e-01,...
    1.218498444843265e+00,...
    1.730095305997005e+00,...
    2.403139445479289e-01,...
    8.526541262207157e-01,...
    1.839467231356713e+00,...
    2.403139445479289e-01,...
    1.350211172437031e+00,...
    9.139950953929686e-01,...
    1.008985635023947e+00]);

H(m-14:m,m-14:m)=rot90(H(1:15,1:15),2);

H=H*h;
HI = inv(H);

a1 = 6/7; a2 = -15/56; a3 = 5/63; a4 = -1/56; a5 = 1/385; a6 = -1/5544;

D1=a6*(diag(ones(m-6,1),6)-diag(ones(m-6,1),-6))+...
    a5*(diag(ones(m-5,1),5)-diag(ones(m-5,1),-5))+...
    a4*(diag(ones(m-4,1),4)-diag(ones(m-4,1),-4))+...
    a3*(diag(ones(m-3,1),3)-diag(ones(m-3,1),-3))+...
    a2*(diag(ones(m-2,1),2)-diag(ones(m-2,1),-2))+...
    a1*(diag(ones(m-1,1),1)-diag(ones(m-1,1),-1));


D1(1,     1)=    -1.735744761135539e+00;
D1(1  ,   2 )=   2.307154400697361e+00;
D1(1,     3)=    -8.049541548921420e-02;
D1(1 ,    4)=    -5.054323197649326e-01;
D1(1 ,   5 )=   -2.276940626773821e-01;
D1(1 ,    6 )=   -4.332046513157440e-03;
D1(1 ,    7 )=   1.922516500571120e-01;
D1(1 ,    8 )=   1.841426475531248e-01;
D1(1 ,    9 )=   4.299540500353117e-02;
D1(1 ,   10 )=   -8.425321510779284e-02;
D1(1 ,   11 )=   -1.783432239793461e-01;
D1(1 ,   12 )=   -9.098340284317097e-02;
D1(1 ,   13 )=   2.989603414709966e-01;
D1(1 ,   14 )=   -1.262666707971657e-01;
D1(1 ,   15 )=   8.040673525575449e-03;
D1(2 ,    1 )=   -4.259233932680148e-01;
D1(2 ,    3  )=  1.077970170802666e-01;
D1(2 ,    4  )=  3.305658398446859e-01;
D1(2 ,    5 )=   1.023936809322525e-01;
D1(2 ,    6  )=  1.004615991085861e-02;
D1(2 ,    7  )=  -9.993439452786020e-02;
D1(2 ,    8  )=  -8.071677622127596e-02;
D1(2 ,    9  )=  -2.794842118590027e-02;
D1(2 ,   10  )=  4.202296437592187e-02;
D1(2 ,   11 )=   8.084986053360355e-02;
D1(2 ,   12 )=   4.970793893896007e-02;
D1(2 ,   13 )=   -1.480838886705940e-01;
D1(2  ,  14  )=  6.449962737470759e-02;
D1(2  ,  15 )=   -5.276215117611654e-03;
D1(3  ,   1  )=  9.648866897889277e-02;
D1(3  ,   2  )=  -6.999339806925928e-01;
D1(3  ,   4  )=  6.447106132171313e-01;
D1(3  ,   5  )=  1.610341017868972e-01;
D1(3  ,   6  )=  -2.926101271261580e-01;
D1(3  ,   7  )=  2.578035970538076e-02;
D1(3  ,   8  )=  -3.649479605631123e-02;
D1(3  ,   9  )=  2.039707558852247e-01;
D1(3  ,  10  )=  -2.443704472345096e-02;
D1(3  ,  11  )=  -4.983260374057472e-02;
D1(3  ,  12  )=  -1.982796172639611e-01;
D1(3  ,  13  )=  2.659121147554237e-01;
D1(3  ,  14  )=  -1.104398155224152e-01;
D1(3  ,  15  )=  1.413137079651327e-02;
D1(4  ,   1  )=  9.507993640872410e-02;
D1(4  ,   2  )=  -3.368441515472113e-01;
D1(4  ,   3  )=  -1.011778679849722e-01;
D1(4  ,   5  )=  1.664700006607916e-01;
D1(4  ,   6  )=  8.420689674052668e-02;
D1(4  ,   7  )=  1.234483107384551e-01;
D1(4  ,   8  )=  3.371285340885326e-02;
D1(4  ,   9  )=  -2.009043511216467e-02;
D1(4  ,  10  )=  -4.311340633236862e-02;
D1(4  ,  11  )=  -3.824719371220764e-02;
D1(4  ,  12  )=  -2.775022659375772e-03;
D1(4  ,  13  )=  8.248678009703486e-02;
D1(4  ,  14  )=  -5.169933586327551e-02;
D1(4  ,  15  )=  8.542635157189999e-03;
D1(5  ,   1  )=  5.267797175762052e-02;
D1(5  ,   2  )=  -1.283203869171410e-01;
D1(5   ,  3  )=  -3.108064515385221e-02;
D1(5  ,   4  )=  -2.047328098663259e-01;
D1(5  ,   6  )=  1.528377576445926e-02;
D1(5  ,   7  )=  1.842140733907868e-01;
D1(5  ,   8  )=  1.776491663357457e-01;
D1(5  ,   9  )=  4.220032524144508e-02;
D1(5  ,  10  )=  -3.611067700750285e-02;
D1(5  ,  11  )=  -1.090730089234403e-01;
D1(5  ,  12  )=  -5.398262792345335e-02;
D1(5  ,  13  )=  1.325672256468355e-01;
D1(5  ,  14  )=  -3.865175590539601e-02;
D1(5  ,  15  )=  -2.640626439781275e-03;
D1(6  ,   1 )=   5.192760351242898e-03;
D1(6  ,   2  )=  -6.523045708997477e-02;
D1(6  ,   3  )=  2.926101271261579e-01;
D1(6  ,   4  )=  -5.365707057867646e-01;
D1(6  ,   5  )=  -7.918783828195401e-02;
D1(6  ,   7  )=  1.455199260610844e-01;
D1(6  ,   8  )=  4.392627009892030e-01;
D1(6  ,   9  )=  -2.877704902282471e-01;
D1(6  ,  10  )=  6.246304940524665e-02;
D1(6  ,  11  )=  -8.757523368284212e-03;
D1(6  ,  12  )=  2.373907173856874e-01;
D1(6  ,  13  )=  -3.887671913712637e-01;
D1(6  ,  14  )=  2.338336212010801e-01;
D1(6  ,  15  )=  -4.998869639321427e-02;
D1( 7 ,    1 )=   -4.544951340626429e-02;
D1( 7 ,    2 )=   1.279732835762292e-01;
D1( 7 ,    3 )=   -5.084438112238577e-03;
D1( 7 ,    4 )=   -1.551380899582667e-01;
D1( 7 ,    5  )=  -1.882368423100825e-01;
D1( 7 ,    6  )=  -2.869964060279157e-02;
D1( 7  ,   8  )=  1.501531969717210e-01;
D1( 7 ,    9  )=  4.398683557974375e-02;
D1( 7  ,  10  )=  9.034991722085359e-02;
D1( 7 ,   11  )=  6.874263504449905e-02;
D1( 7 ,   12  )=  1.191424140454980e-02;
D1( 7 ,   13  )=  -1.387149470006835e-01;
D1( 7,    14  )=  8.081158783540948e-02;
D1( 7 ,   15  )=  -1.260822624267884e-02;
D1( 8 ,    1  )=  -3.065974203065793e-02;
D1( 8  ,   2  )=  7.279861004189266e-02;
D1( 8  ,   3  )=  5.069205358435631e-03;
D1( 8  ,   4  )=  -2.983896481040582e-02;
D1( 8  ,   5  )=  -1.278497691554543e-01;
D1( 8  ,   6  )=  -6.101453024095747e-02;
D1( 8  ,   7  )=  -1.057522301598588e-01;
D1( 8  ,   9  )=  3.518098475279482e-02;
D1( 8  ,  10 )=   1.267951379582364e-01;
D1( 8  ,  11  )=  1.709119775233402e-01;
D1( 8  ,  12 )=   4.871305182138817e-02;
D1( 8  ,  13 )=   -1.381315738947210e-01;
D1( 8  ,  14 )=   2.531339542383916e-02;
D1( 8  ,  15 )=   8.464447412128264e-03;
D1( 9  ,   1  )=  -5.153795872455654e-02;
D1( 9  ,   2  )=  1.814711596347261e-01;
D1( 9  ,   3  )=  -2.039707558852247e-01;
D1( 9  ,   4 )=   1.280172927036423e-01;
D1( 9  ,   5 )=   -2.186470530689342e-01;
D1( 9  ,   6 )=   2.877704902282471e-01;
D1( 9  ,   7 )=   -2.230327950727989e-01;
D1( 9  ,   8 )=   -2.532789210200124e-01;
D1( 9  ,  10 )=   -2.487018120584765e-02;
D1( 9  ,  11 )=   5.656432879788119e-01;
D1( 9  ,  12 )=   -3.179024054582303e-01;
D1( 9  ,  13 )=   3.080272260545759e-01;
D1( 9  ,  14 )=   -2.449533574878925e-01;
D1( 9  ,  15 )=   6.726397132349360e-02;
D1(10  ,   1 )=   2.846411764336633e-02;
D1(10  ,   2 )=   -7.690296670899374e-02;
D1(10  ,   3 )=   6.887391299700927e-03;
D1(10  ,   4 )=   7.742794040563218e-02;
D1(10  ,   5 )=   5.273142829480713e-02;
D1(10  ,   6 )=   -1.760472544430208e-02;
D1(10  ,   7 )=   -1.291159336943500e-01;
D1(10  ,   8  )=  -2.572762697778848e-01;
D1(10  ,   9  )=  7.009467453925059e-03;
D1(10  ,  11  )=  1.983435406810274e-01;
D1(10  ,  12  )=  -3.844771590180474e-02;
D1(10  ,  13  )=  2.943472385247488e-01;
D1(10  ,  14  )=  -1.771193581329739e-01;
D1(10  ,  15  )=  3.146739090039766e-02;
D1(10  ,  16  )=  -2.115455432962849e-04;
D1(11  ,   1  )=  2.792856995883971e-02;
D1(11  ,   2  )=  -6.858300759838749e-02;
D1(11  ,   3  )=  6.510292419375574e-03;
D1(11  ,   4  )=  3.183946691406334e-02;
D1(11  ,   5  )=  7.382985925399284e-02;
D1(11  ,   6  )=  1.144111158506918e-03;
D1(11  ,   7  )=  -4.553644254612264e-02;
D1(11  ,   8  )=  -1.607498111470597e-01;
D1(11  ,   9  )=  -7.389746738841917e-02;
D1(11  ,  10  )=  -9.193881548309495e-02;
D1(11  ,  12  )=  4.858403641270599e-02;
D1(11  ,  13  )=  2.953396945760745e-01;
D1(11  ,  14  )=  -1.897059300798801e-02;
D1(11  ,  15  )=  -2.681387584075025e-02;
D1(11  ,  16  )=  1.412040700223218e-03;
D1(11  ,  17  )=  -9.805838195994568e-05;
D1(12  ,   1  )=  1.090604649488924e-01;
D1(12  ,   2  )=  -3.227573129195589e-01;
D1(12  ,   3  )=  1.982796172639611e-01;
D1(12  ,   4  )=  1.768258806049677e-02;
D1(12  ,   5  )=  2.796931645632894e-01;
D1(12  ,   6  )=  -2.373907173856874e-01;
D1(12  ,   7  )=  -6.041049615427440e-02;
D1(12  ,   8  )=  -3.507005074362792e-01;
D1(12  ,   9  )=  3.179024054582303e-01;
D1(12  ,  10  )=  1.364157359619944e-01;
D1(12  ,  11  )=  -3.718833008893085e-01;
D1(12  ,  13  )=  -1.473587211078091e-01;
D1(12  ,  14  )=  5.107219886498376e-01;
D1(12  ,  15  )=  -1.500513994956780e-02;
D1(12  ,  16  )=  -7.430756001585824e-02;
D1(12  ,  17  )=  1.080837236594302e-02;
D1(12  ,  18  )=  -7.505814143015984e-04;
D1(13  ,   1   )= -6.378169035524298e-02;
D1(13  ,   2   )= 1.711336527390002e-01;
D1(13  ,   3  )=  -4.732769992164881e-02;
D1(13  ,   4  )=  -9.354938176872968e-02;
D1(13  ,   5  )=  -1.222478620034153e-01;
D1(13  ,   6  )=  6.919375218960795e-02;
D1(13  ,   7  )=  1.251833421669687e-01;
D1(13  ,   8  )=  1.769951193441043e-01;
D1(13  ,   9  )=  -5.482345223653074e-02;
D1(13  ,  10  )=  -1.858793591648383e-01;
D1(13  ,  11  )=  -4.023575729350775e-01;
D1(13  ,  12  )=  2.622727189335804e-02;
D1(13  ,  14  )=  4.511356466215911e-01;
D1(13  ,  15  )=  -9.724618589400275e-02;
D1(13  ,  16  )=  5.877975311212341e-02;
D1(13  ,  17  )=  -1.322544445022777e-02;
D1(13  ,  18  )=  1.923701010942221e-03;
D1(13  ,  19  )=  -1.335903479820987e-04;
D1(14   ,  1  )=  3.979504551511218e-02;
D1(14  ,   2  )=  -1.101140274631495e-01;
D1(14  ,   3  )=  2.903760407152547e-02;
D1(14  ,   4  )=  8.661624791356937e-02;
D1(14  ,   5  )=  5.265409651744463e-02;
D1(14  ,   6  )=  -6.148116129069339e-02;
D1(14  ,   7  )=  -1.077344885099467e-01;
D1(14  ,   8  )=  -4.791555974685060e-02;
D1(14  ,   9  )=  6.440483965930399e-02;
D1(14  ,  10  )=  1.652323434850750e-01;
D1(14  ,  11  )=  3.817939983867793e-02;
D1(14  ,  12  )=  -1.342825757801646e-01;
D1(14  ,  13  )=  -6.664460164211102e-01;
D1(14  ,  15  )=  8.751759953804069e-01;
D1(14  ,  16  )=  -2.930619039503475e-01;
D1(14  ,  17  )=  8.683315672602888e-02;
D1(14  ,  18  )=  -1.953746026335650e-02;
D1(14  ,  19  )=  2.841812401942763e-03;
D1(14  ,  20  )=  -1.973480834682474e-04;
D1(15  ,   1  )=  -2.295575530984041e-03;
D1(15  ,   2  )=  8.159561455701054e-03;
D1(15  ,   3  )=  -3.365722305748102e-03;
D1(15  ,   4  )=  -1.296477975670818e-02;
D1(15  ,   5  )=  3.258582992375954e-03;
D1(15  ,   6  )=  1.190599786168119e-02;
D1(15  ,   7  )=  1.522628621820923e-02;
D1(15  ,   8  )=  -1.451388426876270e-02;
D1(15  ,   9  )=  -1.602051576712872e-02;
D1(15  ,  10  )=  -2.659185598017704e-02;
D1(15  ,  11  )=  4.888399224192811e-02;
D1(15  ,  12  )=  3.573831226733737e-03;
D1(15  ,  13  )=  1.301335540499094e-01;
D1(15  ,  14  )=  -7.927829095053140e-01;
D1(15  ,  16  )=  8.495094750506673e-01;
D1(15  ,  17  )=  -2.654717109533336e-01;
D1(15  ,  18  )=  7.865828472691364e-02;
D1(15  ,  19  )=  -1.769811406355557e-02;
D1(15  ,  20  )=  2.574271136517173e-03;
D1(15  ,  21  )=  -1.787688289248037e-04;


D1(m-14:m,m-20:m)=rot90( -D1(1:15,1:21),2);
D1=D1/h;


% Wide 2nd derivative

D2=D1*D1;

e_1 = zeros(m,1);
e_1(1)= 1;
e_m = zeros(m,1);
e_m(end)= 1;

S_1 = (e_1'*D1)';
S_m = (e_m'*D1);

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