annotate +sbp/ordinary10.m @ 87:0a29a60e0b21

In Curve: Rearranged for speed. arc_length_fun is now a property of Curve. If it is not supplied, it is computed via the derivative and spline fitting. Switching to the arc length parameterization is much faster now. The new stuff can be tested with testArcLength.m (which should be deleted after that).
author Martin Almquist <martin.almquist@it.uu.se>
date Sun, 29 Nov 2015 22:23:09 +0100
parents 32b39dc44474
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
29
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
1 function [H, HI, D1, D2, e_1, e_m, M,Q S_1, S_m] = ordinary10(m,h)
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
3 H_U = [0.5261271563e10 / 0.18289152000e11 0 0 0 0 0 0 0 0 0 0; 0 0.2881040311e10 / 0.1828915200e10 0 0 0 0 0 0 0 0 0; 0 0 0.52175551e8 / 0.406425600e9 0 0 0 0 0 0 0 0; 0 0 0 0.11662993e8 / 0.6096384e7 0 0 0 0 0 0 0; 0 0 0 0 0.50124587e8 / 0.87091200e8 0 0 0 0 0 0; 0 0 0 0 0 0.50124587e8 / 0.72576000e8 0 0 0 0 0; 0 0 0 0 0 0 0.148333439e9 / 0.87091200e8 0 0 0 0; 0 0 0 0 0 0 0 0.63867949e8 / 0.152409600e9 0 0 0; 0 0 0 0 0 0 0 0 0.20608675e8 / 0.16257024e8 0 0; 0 0 0 0 0 0 0 0 0 0.1704508063e10 / 0.1828915200e10 0; 0 0 0 0 0 0 0 0 0 0 0.18425967263e11 / 0.18289152000e11;];
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
4
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
5 H=eye(m);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6 H(1:11,1:11)=H_U;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
7 H(m-10:m,m-10:m)=flipud( fliplr(H_U(1:11,1:11) ) );
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8 H=H*h;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9 HI=inv(H);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
10
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
11 Q=-(-1/1260*diag(ones(m-5,1),5)+5/504*diag(ones(m-4,1),4)-5/84*diag(ones(m-3,1),3)+5/21*diag(ones(m-2,1),2)-5/6*diag(ones(m-1,1),1)+5/6*diag(ones(m-1,1),-1)-5/21*diag(ones(m-2,1),-2)+5/84*diag(ones(m-3,1),-3)-5/504*diag(ones(m-4,1),-4)+1/1260*diag(ones(m-5,1),-5));
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
12
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
13 Q_U = [-0.1e1 / 0.2e1 0.2300876759589119e16 / 0.3395198177280000e16 -0.99808615498093e14 / 0.2263465451520000e16 -0.34957747037683e14 / 0.212199886080000e15 -0.709586095717e12 / 0.13473008640000e14 0.325330433051e12 / 0.6218311680000e13 0.27953548723573e14 / 0.485028311040000e15 0.2690678501e10 / 0.412439040000e12 -0.2397491025029e13 / 0.70733295360000e14 -0.9959492094287e13 / 0.1131732725760000e16 0.5242772857661e13 / 0.522338181120000e15; -0.2300876759589119e16 / 0.3395198177280000e16 0 0.3103439505511e13 / 0.16643128320000e14 0.2700334568377e13 / 0.5052378240000e13 0.50587599589937e14 / 0.242514155520000e15 -0.5570893587157e13 / 0.40419025920000e14 -0.1496329934863e13 / 0.8083805184000e13 -0.322512443237e12 / 0.12482346240000e14 0.2275340833763e13 / 0.23096586240000e14 0.22922115021893e14 / 0.848799544320000e15 -0.143e3 / 0.5000e4; 0.99808615498093e14 / 0.2263465451520000e16 -0.3103439505511e13 / 0.16643128320000e14 0 0.15053664233879e14 / 0.40419025920000e14 -0.9306441440587e13 / 0.32335220736000e14 -0.945459729233e12 / 0.13473008640000e14 0.956829267413e12 / 0.5774146560000e13 0.446866085903e12 / 0.56586636288000e14 -0.41109372242993e14 / 0.754488483840000e15 0.1e1 / 0.500e3 0.17e2 / 0.2500e4; 0.34957747037683e14 / 0.212199886080000e15 -0.2700334568377e13 / 0.5052378240000e13 -0.15053664233879e14 / 0.40419025920000e14 0 0.3899174751943e13 / 0.10104756480000e14 0.4691717443831e13 / 0.10104756480000e14 -0.58571891887e11 / 0.396264960000e12 0.100791910589e12 / 0.1040195520000e13 -0.425149181e9 / 0.29719872000e11 -0.2376515922259e13 / 0.30314269440000e14 0.36894656431e11 / 0.1036385280000e13; 0.709586095717e12 / 0.13473008640000e14 -0.50587599589937e14 / 0.242514155520000e15 0.9306441440587e13 / 0.32335220736000e14 -0.3899174751943e13 / 0.10104756480000e14 0 -0.4552305973e10 / 0.444165120000e12 0.4984940784247e13 / 0.11548293120000e14 -0.19410791e8 / 0.146764800e9 -0.2912773695913e13 / 0.40419025920000e14 0.127067639161e12 / 0.3233522073600e13 -0.89277540287e11 / 0.37309870080000e14; -0.325330433051e12 / 0.6218311680000e13 0.5570893587157e13 / 0.40419025920000e14 0.945459729233e12 / 0.13473008640000e14 -0.4691717443831e13 / 0.10104756480000e14 0.4552305973e10 / 0.444165120000e12 0 0.31722122083e11 / 0.84913920000e11 -0.887187251021e12 / 0.10104756480000e14 -0.1661755478749e13 / 0.26946017280000e14 0.1505713246249e13 / 0.13473008640000e14 -0.38859042469e11 / 0.1036385280000e13; -0.27953548723573e14 / 0.485028311040000e15 0.1496329934863e13 / 0.8083805184000e13 -0.956829267413e12 / 0.5774146560000e13 0.58571891887e11 / 0.396264960000e12 -0.4984940784247e13 / 0.11548293120000e14 -0.31722122083e11 / 0.84913920000e11 0 0.9357094407023e13 / 0.20209512960000e14 0.52602356173249e14 / 0.161676103680000e15 -0.1435252677707e13 / 0.17322439680000e14 -0.33048158431e11 / 0.3109155840000e13; -0.2690678501e10 / 0.412439040000e12 0.322512443237e12 / 0.12482346240000e14 -0.446866085903e12 / 0.56586636288000e14 -0.100791910589e12 / 0.1040195520000e13 0.19410791e8 / 0.146764800e9 0.887187251021e12 / 0.10104756480000e14 -0.9357094407023e13 / 0.20209512960000e14 0 0.70089734285659e14 / 0.141466590720000e15 -0.105938137621e12 / 0.471555302400e12 0.4358988450443e13 / 0.65292272640000e14; 0.2397491025029e13 / 0.70733295360000e14 -0.2275340833763e13 / 0.23096586240000e14 0.41109372242993e14 / 0.754488483840000e15 0.425149181e9 / 0.29719872000e11 0.2912773695913e13 / 0.40419025920000e14 0.1661755478749e13 / 0.26946017280000e14 -0.52602356173249e14 / 0.161676103680000e15 -0.70089734285659e14 / 0.141466590720000e15 0 0.314274398580227e15 / 0.377244241920000e15 -0.97822819709e11 / 0.487710720000e12; 0.9959492094287e13 / 0.1131732725760000e16 -0.22922115021893e14 / 0.848799544320000e15 -0.1e1 / 0.500e3 0.2376515922259e13 / 0.30314269440000e14 -0.127067639161e12 / 0.3233522073600e13 -0.1505713246249e13 / 0.13473008640000e14 0.1435252677707e13 / 0.17322439680000e14 0.105938137621e12 / 0.471555302400e12 -0.314274398580227e15 / 0.377244241920000e15 0 0.7519148725913e13 / 0.9327467520000e13; -0.5242772857661e13 / 0.522338181120000e15 0.143e3 / 0.5000e4 -0.17e2 / 0.2500e4 -0.36894656431e11 / 0.1036385280000e13 0.89277540287e11 / 0.37309870080000e14 0.38859042469e11 / 0.1036385280000e13 0.33048158431e11 / 0.3109155840000e13 -0.4358988450443e13 / 0.65292272640000e14 0.97822819709e11 / 0.487710720000e12 -0.7519148725913e13 / 0.9327467520000e13 0;];
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
14
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
15 Q(1:11,1:11)=Q_U;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
16 Q(m-10:m,m-10:m)=flipud( fliplr(-Q_U ) );
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
17
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
18 D1=HI*Q;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
19
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
20
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
21 s18=-1.000000000000000; s19=0.195000000000000; % alpha 0.0605
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
22 s18= -0.475000000000000; s19=0.110000000000000; % alpha 0.0350
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
23 %s18=0;s19=0;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
24 DS=zeros(m,m);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
25 DS(1,1:9)=[0.49e2 / 0.20e2 - s18 - (7 * s19) -0.6e1 + 0.7e1 * s18 + (48 * s19) 0.15e2 / 0.2e1 - 0.21e2 * s18 - (140 * s19) -0.20e2 / 0.3e1 + 0.35e2 * s18 + (224 * s19) 0.15e2 / 0.4e1 - 0.35e2 * s18 - (210 * s19) -0.6e1 / 0.5e1 + 0.21e2 * s18 + (112 * s19) 0.1e1 / 0.6e1 - 0.7e1 * s18 - (28 * s19) s18 s19;];
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
26 DS(m,m-8:m)=fliplr(DS(1,1:9));
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
27 DS=DS/h;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
28
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
29 M_U = [0.12056593789671863908e1 -0.13378814169347239658e1 0.36847309286546532061e-2 0.15698288365600946515e0 -0.37472461482539197952e-2 -0.62491712449361657064e-2 -0.29164045872729581661e-1 0.54848184117832929161e-3 0.13613461413384884448e-1 -0.25059220258337808220e-2 -0.94113457993630916498e-3; -0.13378814169347239658e1 0.21749807117105597139e1 -0.12369059547124894597e0 -0.83712574037924152603e0 0.50065127254670973258e-1 0.81045853127317536361e-2 0.97405846039248226536e-1 -0.68942461520402214720e-3 -0.41326971493379188475e-1 0.75778529605774119402e-2 0.25800256160095691057e-2; 0.36847309286546532061e-2 -0.12369059547124894597e0 0.18361596652499065332e0 0.48289690013342693109e-1 -0.19719621435164680412e0 0.11406859029505842791e0 -0.29646295985488126964e-1 -0.16038463172861201306e-2 0.32879841528337653050e-2 -0.93242311589807387463e-3 0.12241332668787820533e-3; 0.15698288365600946515e0 -0.83712574037924152603e0 0.48289690013342693109e-1 0.12886524606662484673e1 -0.14403037739488789185e0 -0.44846291607489015475e0 -0.10598334599408054277e0 -0.15873275740355918053e-1 0.73988493386459608166e-1 -0.12508848749152899785e-1 -0.39290233894513005339e-2; -0.37472461482539197952e-2 0.50065127254670973258e-1 -0.19719621435164680412e0 -0.14403037739488789185e0 0.51482665719685186210e0 0.51199577887125103015e-1 -0.36233561810883077365e0 0.91356850268746392169e-1 0.24195916108052419451e-2 -0.18564214413731389338e-2 -0.70192677320704413827e-3; -0.62491712449361657064e-2 0.81045853127317536361e-2 0.11406859029505842791e0 -0.44846291607489015475e0 0.51199577887125103015e-1 0.68636003380365860083e0 -0.28358848290867614908e0 -0.13836006478253396528e0 0.76158070663111995297e-2 0.11447010307180005164e-1 -0.21349696610286552676e-2; -0.29164045872729581661e-1 0.97405846039248226536e-1 -0.29646295985488126964e-1 -0.10598334599408054277e0 -0.36233561810883077365e0 -0.28358848290867614908e0 0.15216081480839085990e1 -0.42653865162216293237e0 -0.42047484981879143123e0 0.19813359263872926304e-1 0.19221397241190103344e-1; 0.54848184117832929161e-3 -0.68942461520402214720e-3 -0.16038463172861201306e-2 -0.15873275740355918053e-1 0.91356850268746392169e-1 -0.13836006478253396528e0 -0.42653865162216293237e0 0.10656733504627815335e1 -0.66921872668484232217e0 0.12022033144141336599e0 -0.30157881394591483631e-1; 0.13613461413384884448e-1 -0.41326971493379188475e-1 0.32879841528337653050e-2 0.73988493386459608166e-1 0.24195916108052419451e-2 0.76158070663111995297e-2 -0.42047484981879143123e0 -0.66921872668484232217e0 0.24064247712949611684e1 -0.15150200315922263367e1 0.17373015320416595052e0; -0.25059220258337808220e-2 0.75778529605774119402e-2 -0.93242311589807387463e-3 -0.12508848749152899785e-1 -0.18564214413731389338e-2 0.11447010307180005164e-1 0.19813359263872926304e-1 0.12022033144141336599e0 -0.15150200315922263367e1 0.27682502485427255096e1 -0.15975407111468405444e1; -0.94113457993630916498e-3 0.25800256160095691057e-2 0.12241332668787820533e-3 -0.39290233894513005339e-2 -0.70192677320704413827e-3 -0.21349696610286552676e-2 0.19221397241190103344e-1 -0.30157881394591483631e-1 0.17373015320416595052e0 -0.15975407111468405444e1 0.29033627686681129471e1;];
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
30
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
31
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
32 M=-(1/3150)*diag(ones(m-5,1),5)+(5/1008)*diag(ones(m-4,1),4)-(5/126)*diag(ones(m-3,1),3)+(5/21)*diag(ones(m-2,1),2)-(5/3)*diag(ones(m-1,1),1)...
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
33 -(1/3150)*diag(ones(m-5,1),-5)+(5/1008)*diag(ones(m-4,1),-4)-(5/126)*diag(ones(m-3,1),-3)+(5/21)*diag(ones(m-2,1),-2)-(5/3)*diag(ones(m-1,1),-1)...
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
34 +(5269/1800)*diag(ones(m,1),0);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
35
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
36 M(1:11,1:11)=M_U;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
37 M(m-10:m,m-10:m)=flipud( fliplr(M_U(1:11,1:11) ) );
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
38
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
39 D2=HI*(-M/h+DS);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
40
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
41 e_1 = zeros(m,1);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
42 e_1(1)= 1;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
43 e_m = zeros(m,1);
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
44 e_m(end)= 1;
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
45 S_1 = -DS(1,:)';
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
46 S_m = DS(end,:)';
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
47
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
48 Q = H*D1-(-e_1*e_1' + e_m*e_m');
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
49 M = -(H*D2-(-e_1*S_1' + e_m*S_m'));
32b39dc44474 Removed repository inside +sbp to make it part of the root repo.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
50 end