Mercurial > matrix-functions
view 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 source
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;