Mercurial > matrix-functions
diff mftoolbox/cosm_pade.m @ 0:8f23314345f4 draft
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 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mftoolbox/cosm_pade.m Wed May 06 14:56:53 2015 +0200 @@ -0,0 +1,122 @@ +function C = cosm_pade(A,m,sq) +%COSM_PADE Evaluate Pade approximation to the matrix cosine. +% C = COSM_PADE(A,M,SQ) approximates the matrix cosine using the +% Mth order diagonal Pade approximation. +% If SQ = 1 (default) then C is an approximation to cos(sqrt(A)); +% otherwise C is an approximation to cos(A). + +if nargin < 3 + sq = 1; +end +if sq == 1 + A2 = A; +else + A2 = A^2; +end + +n = length(A2); +I = eye(n); +if m == 2 + X2 = A2; + P = I - (5/12)*X2; + Q = I + (1/12)*X2; +elseif m == 4 + X2 = A2; X4 = X2^2; + P = I - (115/252)*X2 + (313/15120)*X4; + Q = I + (11/252)*X2 + (13/15120)*X4; +elseif m == 6 + X2 = A2; X4 = X2^2; X6 = X4*X2; + P = I - (3665/7788)*X2 + (711/25960)*X4 - (2923/7850304)*X6; + Q = I + (229/7788)*X2 + (1/2360)*X4 + (127/39251520)*X6; +elseif m == 8 + X2 = A2; X4 = X2^2; X6 = X4*X2; X8 = X6*X2; + P = I - (260735/545628)*X2 + (4375409/141863280)*X4 ... + - (7696415/13108167072)*X6 + (80737373/23594700729600)*X8 ; + Q = I + (12079/545628)*X2 + (34709/141863280)*X4 ... + + (109247/65540835360)*X6 + (11321/1814976979200)*X8 ; +elseif m == 12 + X1 = A2; X2 = X1*X1; X3 = X2*X1; + p = [1,-220574348151635/454605030049116,20837207639809/606140040065488,... + -199961484798769/241849875986129712,38062401688454831/... + 4440363723105341512320,-116112688080827/2894459315802000393216,... + 151259208063389819/2133505961677654489839513600]; + q = [1,6728166872923/454605030049116,66817219029/606140040065488,... + 650617920073/1209249379930648560,8225608067111/4440363723105341512320,... + 2848116281867/651253346055450088473600,12170851069679/... + 2133505961677654489839513600]; + P = X3*((p(7)*eye(n))*X3+(p(4)*eye(n)+p(5)*X1+p(6)*X2)*eye(n))+... + (p(1)*eye(n)+p(2)*X1+p(3)*X2); + Q = X3*((q(7)*eye(n))*X3+(q(4)*eye(n)+q(5)*X1+q(6)*X2)*eye(n))+... + (q(1)*eye(n)+q(2)*X1+q(3)*X2); +elseif m == 16 + X1 = A2; X2 = X1*X1; X3 = X2*X1; X4 = X2*X2; + p = [1,-1126682407530029115789472765/2304577612359442026681336372,... + 145053661043845297596963732421/4009965045505429126425525287280,... + -1534672316720770887322573595/1603986018202171650570210114912,... + 718202654899849477670594159641/60630671488042088391553942343673600,... + -128936233968950140829066659951/1673406533069961639606888808685391360,... + 6524116556754642812271854422129/23929713422900451446378509964201096448000,... + -382586638331055978467487427009/... + 763836452458982410168402038057298998620160,... + 88555612088268453352055067469523/... + 233733954452448617511531023645533493577768960000]; + q = [1,25606398649691897551195421/2304577612359442026681336372,... + 22668270274336502918805611/364542276864129920584138662480,... + 1853378279158412863783499/8019930091010858252851050574560,... + 38226389122327179481602241/60630671488042088391553942343673600,... + 995615371594253927197913/760639333213618927094040367584268800,... + 225870994754204367988837/110275177064057379937228156517055744000,... + 42889724495628101076622829/19095911311474560254210050951432474965504000,... + 2603898999593850290644763/1931685573987178657120091104508541269237760000]; + P = X4*((p(9)*eye(n))*X4+(p(5)*eye(n)+p(6)*X1+p(7)*X2+p(8)*X3)*eye(n))+... + (p(1)*eye(n)+p(2)*X1+p(3)*X2+p(4)*X3); + Q = X4*((q(9)*eye(n))*X4+(q(5)*eye(n)+q(6)*X1+q(7)*X2+q(8)*X3)*eye(n))+... + (q(1)*eye(n)+q(2)*X1+q(3)*X2+q(4)*X3); +elseif m == 20 + X1 = A2; X2 = X1*X1; X3 = X2*X1; X4 = X2*X2; X5 = X4*X1; + p = [1,-18866133841442352341137832915472113127673/... + 38415527280635118612047973206722428679860,... + 917980006162069077942240197016800349995791/... + 24637158162647322736526766816577984260016880,... + -4028339250935885155796261896908967142863591/... + 3880352410616953331002965773611032520952658600,... + 3925400573997340625949450726927185904756763/... + 279385373564420639832213535699994341508591419200,... + -804035081520215224783821741744290679884325097/... + 7621632990837395054622785253895845636354373915776000,... + 19795406323827219175300218252334434555489703/... + 42100448901768467920773480450091337800814636868096000,... + -118523567829079039162888326742509818128969627/... + 92831489828399471765305524392451399850796274294151680000,... + 6151694105279089780298999575203793198700323571/... + 2954269332298984789459083008265373348851740633137083064320000,... + -438673281197605688527510681818034658057709668453/... + 232382825678638143538851469430154267620677918202562953839411200000,... + 31699084606166905465868332652040368902350407479/... + 42944346185412328925979751550692508656301279283833633869523189760000]; + q = [1,341629798875206964886153687889101212257/... + 38415527280635118612047973206722428679860,... + 981038224413663993784862242489461225499/... + 24637158162647322736526766816577984260016880,... + 461441299765418864926911910257258436499/... + 3880352410616953331002965773611032520952658600,... + 73764947345500690357380325430051300659/... + 279385373564420639832213535699994341508591419200,... + 3496016725011957790816142668159659762953/... + 7621632990837395054622785253895845636354373915776000,... + 562876526229442596170390468670872658343/... + 884109426937137826336243089451918093817107374230016000,... + 65309174262483666596220950666851746623/... + 92831489828399471765305524392451399850796274294151680000,... + 1768262649350763278383509302712194678051/... + 2954269332298984789459083008265373348851740633137083064320000,... + 83263779334467686055536878437959026858717/... + 232382825678638143538851469430154267620677918202562953839411200000,... + 24953265550459114615706087077367245444511/... + 214721730927061644629898757753462543281506396419168169347615948800000]; + 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))... + +p(1)*eye(n)+p(2)*X1+p(3)*X2+p(4)*X3+p(5)*X4; + 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))... + +(q(1)*eye(n)+q(2)*X1+q(3)*X2+q(4)*X3+q(5)*X4); +end +C = Q\P;