Mercurial > repos > public > sbplib

comparison +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
comparison
equal deleted inserted replaced
266:bfa130b7abf6 267:f7ac3cd6eeaa
4 BP = 8; 4 BP = 8;
5 if(m<2*BP) 5 if(m<2*BP)
6 error(['Operator requires at least ' num2str(2*BP) ' grid points']); 6 error(['Operator requires at least ' num2str(2*BP) ' grid points']);
7 end 7 end
8 8
9 e_1=zeros(m,1);e_1(1)=1; 9 e_1=sparse(m,1);e_1(1)=1;
10 e_m=zeros(m,1);e_m(m)=1; 10 e_m=sparse(m,1);e_m(m)=1;
11 11
12 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; 12 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;
13 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; 13 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;
14 -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; 14 -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;
15 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; 15 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;
18 -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; 18 -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;
19 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; 19 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;
20 ]; 20 ];
21 21
22 22
23 H=eye(m); 23 H=speye(m);
24 H(1:8,1:8)=H_U; 24 H(1:8,1:8)=H_U;
25 H(m-7:m,m-7:m)=rot90( H_U(1:8,1:8) ,2 ); 25 H(m-7:m,m-7:m)=rot90( H_U(1:8,1:8) ,2 );
26 H=H*h; 26 H=H*h;
27 HI=inv(H); 27 HI=inv(H);
28 28
29 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)); 29 % 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));
30 30 diags = -4:4;
31 stencil = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280];
32 Q = stripeMatrix(stencil, diags, m);
31 33
32 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;]; 34 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;];
33 35
34 Q(1:8,1:8)=Q_U; 36 Q(1:8,1:8)=Q_U;
35 Q(m-7:m,m-7:m)=rot90( -Q_U(1:8,1:8) ,2 ); 37 Q(m-7:m,m-7:m)=rot90( -Q_U(1:8,1:8) ,2 );
36 38
37 D1=HI*Q; 39 D1=H\Q;
38 40
39 41
40 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;]; 42 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;];
41 43
42 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)); 44 % 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));
45 diags = -4:4;
46 left_stencil = [-1/560,8/315,-1/5,8/5];
47 stencil = [left_stencil,-205/72,fliplr(left_stencil)];
48 M = stripeMatrix(stencil, diags, m);
49
43 M(1:8,1:8)=M_U; 50 M(1:8,1:8)=M_U;
44 51
45 M(m-7:m,m-7:m)=rot90( M_U ,2 ); 52 M(m-7:m,m-7:m)=rot90( M_U ,2 );
46 M=M/h; 53 M=M/h;
47 54
48 % 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;]; 55 % 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;];
49 % DS=zeros(m,m); 56 % DS=sparse(m,m);
50 % DS(1,1:8)=DS_U; 57 % DS(1,1:8)=DS_U;
51 % DS(m,m-7:m)=fliplr(DS_U); 58 % DS(m,m-7:m)=fliplr(DS_U);
52 % DS=DS/h; 59 % DS=DS/h;
53 60
54 %D2=HI*(-M+DS); 61 %D2=HI*(-M+DS);
55 62
56 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; 63 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;
57 S_1=zeros(1,m); 64 S_1=sparse(1,m);
58 S_1(1:8)=S_U; 65 S_1(1:8)=S_U;
59 S_m=zeros(1,m); 66 S_m=sparse(1,m);
60 S_m(m-7:m)=fliplr(-S_U); 67 S_m(m-7:m)=fliplr(-S_U);
61 68
62 69
63 D2=HI*(-M - e_1*S_1+e_m*S_m); 70 D2=H\(-M - e_1*S_1+e_m*S_m);
64 S_1 = S_1'; 71 S_1 = S_1';
65 S_m = S_m'; 72 S_m = S_m';
66 end 73 end