Mercurial > matrix-functions

comparison mftoolbox/cosm_pade.m @ 0:8f23314345f4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Create local repository for matrix toolboxes. Step #0 done.
author Antonio Pino Robles <data.script93@gmail.com>
date Wed, 06 May 2015 14:56:53 +0200
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:8f23314345f4
1 function C = cosm_pade(A,m,sq)
2 %COSM_PADE Evaluate Pade approximation to the matrix cosine.
3 % C = COSM_PADE(A,M,SQ) approximates the matrix cosine using the
4 % Mth order diagonal Pade approximation.
5 % If SQ = 1 (default) then C is an approximation to cos(sqrt(A));
6 % otherwise C is an approximation to cos(A).
7
8 if nargin < 3
9 sq = 1;
10 end
11 if sq == 1
12 A2 = A;
13 else
14 A2 = A^2;
15 end
16
17 n = length(A2);
18 I = eye(n);
19 if m == 2
20 X2 = A2;
21 P = I - (5/12)*X2;
22 Q = I + (1/12)*X2;
23 elseif m == 4
24 X2 = A2; X4 = X2^2;
25 P = I - (115/252)*X2 + (313/15120)*X4;
26 Q = I + (11/252)*X2 + (13/15120)*X4;
27 elseif m == 6
28 X2 = A2; X4 = X2^2; X6 = X4*X2;
29 P = I - (3665/7788)*X2 + (711/25960)*X4 - (2923/7850304)*X6;
30 Q = I + (229/7788)*X2 + (1/2360)*X4 + (127/39251520)*X6;
31 elseif m == 8
32 X2 = A2; X4 = X2^2; X6 = X4*X2; X8 = X6*X2;
33 P = I - (260735/545628)*X2 + (4375409/141863280)*X4 ...
34 - (7696415/13108167072)*X6 + (80737373/23594700729600)*X8 ;
35 Q = I + (12079/545628)*X2 + (34709/141863280)*X4 ...
36 + (109247/65540835360)*X6 + (11321/1814976979200)*X8 ;
37 elseif m == 12
38 X1 = A2; X2 = X1*X1; X3 = X2*X1;
39 p = [1,-220574348151635/454605030049116,20837207639809/606140040065488,...
40 -199961484798769/241849875986129712,38062401688454831/...
41 4440363723105341512320,-116112688080827/2894459315802000393216,...
42 151259208063389819/2133505961677654489839513600];
43 q = [1,6728166872923/454605030049116,66817219029/606140040065488,...
44 650617920073/1209249379930648560,8225608067111/4440363723105341512320,...
45 2848116281867/651253346055450088473600,12170851069679/...
46 2133505961677654489839513600];
47 P = X3*((p(7)*eye(n))*X3+(p(4)*eye(n)+p(5)*X1+p(6)*X2)*eye(n))+...
48 (p(1)*eye(n)+p(2)*X1+p(3)*X2);
49 Q = X3*((q(7)*eye(n))*X3+(q(4)*eye(n)+q(5)*X1+q(6)*X2)*eye(n))+...
50 (q(1)*eye(n)+q(2)*X1+q(3)*X2);
51 elseif m == 16
52 X1 = A2; X2 = X1*X1; X3 = X2*X1; X4 = X2*X2;
53 p = [1,-1126682407530029115789472765/2304577612359442026681336372,...
54 145053661043845297596963732421/4009965045505429126425525287280,...
55 -1534672316720770887322573595/1603986018202171650570210114912,...
56 718202654899849477670594159641/60630671488042088391553942343673600,...
57 -128936233968950140829066659951/1673406533069961639606888808685391360,...
58 6524116556754642812271854422129/23929713422900451446378509964201096448000,...
59 -382586638331055978467487427009/...
60 763836452458982410168402038057298998620160,...
61 88555612088268453352055067469523/...
62 233733954452448617511531023645533493577768960000];
63 q = [1,25606398649691897551195421/2304577612359442026681336372,...
64 22668270274336502918805611/364542276864129920584138662480,...
65 1853378279158412863783499/8019930091010858252851050574560,...
66 38226389122327179481602241/60630671488042088391553942343673600,...
67 995615371594253927197913/760639333213618927094040367584268800,...
68 225870994754204367988837/110275177064057379937228156517055744000,...
69 42889724495628101076622829/19095911311474560254210050951432474965504000,...
70 2603898999593850290644763/1931685573987178657120091104508541269237760000];
71 P = X4*((p(9)*eye(n))*X4+(p(5)*eye(n)+p(6)*X1+p(7)*X2+p(8)*X3)*eye(n))+...
72 (p(1)*eye(n)+p(2)*X1+p(3)*X2+p(4)*X3);
73 Q = X4*((q(9)*eye(n))*X4+(q(5)*eye(n)+q(6)*X1+q(7)*X2+q(8)*X3)*eye(n))+...
74 (q(1)*eye(n)+q(2)*X1+q(3)*X2+q(4)*X3);
75 elseif m == 20
76 X1 = A2; X2 = X1*X1; X3 = X2*X1; X4 = X2*X2; X5 = X4*X1;
77 p = [1,-18866133841442352341137832915472113127673/...
78 38415527280635118612047973206722428679860,...
79 917980006162069077942240197016800349995791/...
80 24637158162647322736526766816577984260016880,...
81 -4028339250935885155796261896908967142863591/...
82 3880352410616953331002965773611032520952658600,...
83 3925400573997340625949450726927185904756763/...
84 279385373564420639832213535699994341508591419200,...
85 -804035081520215224783821741744290679884325097/...
86 7621632990837395054622785253895845636354373915776000,...
87 19795406323827219175300218252334434555489703/...
88 42100448901768467920773480450091337800814636868096000,...
89 -118523567829079039162888326742509818128969627/...
90 92831489828399471765305524392451399850796274294151680000,...
91 6151694105279089780298999575203793198700323571/...
92 2954269332298984789459083008265373348851740633137083064320000,...
93 -438673281197605688527510681818034658057709668453/...
94 232382825678638143538851469430154267620677918202562953839411200000,...
95 31699084606166905465868332652040368902350407479/...
96 42944346185412328925979751550692508656301279283833633869523189760000];
97 q = [1,341629798875206964886153687889101212257/...
98 38415527280635118612047973206722428679860,...
99 981038224413663993784862242489461225499/...
100 24637158162647322736526766816577984260016880,...
101 461441299765418864926911910257258436499/...
102 3880352410616953331002965773611032520952658600,...
103 73764947345500690357380325430051300659/...
104 279385373564420639832213535699994341508591419200,...
105 3496016725011957790816142668159659762953/...
106 7621632990837395054622785253895845636354373915776000,...
107 562876526229442596170390468670872658343/...
108 884109426937137826336243089451918093817107374230016000,...
109 65309174262483666596220950666851746623/...
110 92831489828399471765305524392451399850796274294151680000,...
111 1768262649350763278383509302712194678051/...
112 2954269332298984789459083008265373348851740633137083064320000,...
113 83263779334467686055536878437959026858717/...
114 232382825678638143538851469430154267620677918202562953839411200000,...
115 24953265550459114615706087077367245444511/...
116 214721730927061644629898757753462543281506396419168169347615948800000];
117 P = X5*((p(11)*eye(n))*X5+(p(6)*eye(n)+p(7)*X1+p(8)*X2+p(9)*X3+p(10)*X4)*eye(n))...
118 +p(1)*eye(n)+p(2)*X1+p(3)*X2+p(4)*X3+p(5)*X4;
119 Q = X5*((q(11)*eye(n))*X5+(q(6)*eye(n)+q(7)*X1+q(8)*X2+q(9)*X3+q(10)*X4)*eye(n))...
120 +(q(1)*eye(n)+q(2)*X1+q(3)*X2+q(4)*X3+q(5)*X4);
121 end
122 C = Q\P;