www.gusucode.com > matlab矩阵函数工具箱 > matlab矩阵函数工具箱/matlab矩阵函数工具箱/cosmsinm_pade.m
function [C,S] = cosmsinm_pade(A,m)
%COSMSINM_PADE Evaluate Pade approximations to matrix cosine and sine.
% [C,S] = COSMSINM_PADE(A,M) approximates the matrix cosine
% C = cos(A) and the matrix sine S = sin(A) using the [M/M] and
% [M+1/M+1] Pade approximants, respectively, where M is even.
n = length(A);
I = eye(n);
if m == 2
X2 = A^2;
p = I - (5/12)*X2;
q = I + (1/12)*X2;
p2 = A*(-7/60*X2+I);
q2 = I+1/20*X2;
elseif m == 4
X2 = A^2; X4 = X2^2;
p = I - (115/252)*X2 + (313/15120)*X4;
q = I + (11/252)*X2 + (13/15120)*X4;
p2 = A*(551/166320*X4-53/396*X2+I);
q2 = I+13/396*X2+5/11088*X4;
elseif m == 6
X2 = A^2; 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;
p2 = A*(-479249/11511339840*X6+34911/7613320*X4-29593/207636*X2+I);
q2 = I+1671/69212*X2+97/351384*X4+2623/1644477120*X6;
elseif m == 8
X2 = A^2; 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 ;
p2 = A*(4585922449/15605159573203200*X8-269197963/3940696861920*X6...
+ 38518909/7217393520*X4-53272705/360869676*X2+I);
q2 = I+2290747/120289892*X2+1281433/7217393520*X4+560401/...
562956694560*X6+1029037/346781323848960*X8;
elseif m == 10
X2 = A^2; X4 = X2^2; X6 = X4*X2; X8 = X6*X2; X10 = X8*X2;
p = I-213692663231/11226787947026555904*X10+2677576097699/...
425257119205551360*X8-18894084119/25961973089472*X6+...
53471234645/1622623318092*X4-5114526085/10605381164*X2;
q = I+2045322787/280669698675663897600*X10+1469628299/425257119205551360*...
X8+117715523/129809865447360*X6+256513745/1622623318092*X4+...
188164497/10605381164*X2;
p2 = A*(-481959816488503/363275871831577908403200*X10+23704595077729/...
42339845201815607040*X8-2933434889971/33603051747472704*X6+...
3605886663403/617703157122660*X4-109061004303/722459832892*X2+I);
q2 = I+34046903537/2167379498676*X2+1679739379/13726736824948*X4+...
101555058991/168015258737363520*X6+3924840709/2016183104848362240*X8+...
37291724011/11008359752472057830400*X10;
elseif m == 12
X2 = A^2; X4 = X2^2; X6 = X4*X2; X8 = X6*X2; X10 = X8*X2; X12 = X10*X2;
p = I-220574348151635/454605030049116*X2+20837207639809/606140040065488*X4-...
199961484798769/241849875986129712*X6+38062401688454831/...
4440363723105341512320*X8-116112688080827/2894459315802000393216*X10+...
151259208063389819/2133505961677654489839513600*X12;
q = I+6728166872923/454605030049116*X2+66817219029/606140040065488*X4+...
650617920073/1209249379930648560*X6+8225608067111/...
4440363723105341512320*X8+2848116281867/651253346055450088473600*X10+...
12170851069679/2133505961677654489839513600*X12;
p2 = A*(29255060349997508141/7036046883826051377386580480000*X12-...
2614538885979566893/906830082707196962219136000*X10+366447091815846989/...
467616734780258173795200*X8-404844663233524697/...
3987959482775106258000*X6+235300071854130877/37980566502620059600*X4-...
189866343920582047/1238496733781088900*X2+I);
q2 = I+5516592792088701/412832244593696300*X2+10149321432971387/...
113941699507860178800*X4+307158833383841/797591896555021251600*X6+...
234442582333727/202056613793938717072000*X8+5257566411989069/...
2225855657554028907265152000*X10+547584669013343/...
209866390569379868399285760000*X12;
elseif m == 14
X2 = A^2; X4 = X2^2; X6 = X4*X2; X8 = X6*X2; X10 = X8*X2; X12 = X10*X2; ...
X14 = X12*X2;
p = I-9494011832127130075/19482625349840916996*X2+310367939726544641761/...
8767181407428412648200*X4-290442017949426019417/...
322632275793365585453760*X6+422064995289725174621/...
40651666749964063767173760*X8-10143257867238927098753/...
169923967014849786546786316800*X10+569672335133865317479319/...
3379787703925362254415579841152000*X12-42948386132766526786783/...
227121733703784343496726965325414400*X14;
q = I+247300842793328423/19482625349840916996*X2+711404045526343261/...
8767181407428412648200*X4+184363345338626039/...
537720459655609309089600*X6+42648708634036181/...
40651666749964063767173760*X8+2008837845325770881/...
849619835074248932733931584000*X10+1790291961275627717/...
482826814846480322059368548736000*X12+1651330565298296101/...
516185758417691689765288557557760000*X14;
p2 = A*(-91320110661332736729305509/...
9540763134798175069489060319738357760000*X14+...
261319443906679357131557951/25722645706563707295191093999294592000*X12...
-25455044236652128985986613/5785570334359808208545005397952000*X10+...
145834078116038838318251/150274554139215797624545594752*X8-...
264802008995134681176249/2352661207225153389183707200*X6+...
1466350032750183270853069/226863759268139791099857480*X4-...
390816793511920445322905/2520708436312664345553972*X2+I);
q2 = I+9767093068952315200919/840236145437554781851324*X2+...
3067578723707839145789/45372751853627958219971496*X4+...
7023288908523954335543/27223651112176774931982897600*X6+...
2065373697229749176197/2922005219373640509366164342400*X8+...
106107024079408237517/75137277069607898812272797376000*X10+...
1279669772129105715749/659555018117018135774130615366528000*X12+...
70807382924337520289663/48975917425297298690043842974656903168000*X14;
elseif m == 16
X1 = A^2; 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];
p2 = fliplr([8061385294694990486176108590199657/...
477924965253252594838839715519805972908219146240000,...
-4615365296919423054990558658047229/...
177615978658790247602767966416950841324802560000,...
99552990732049352416683777921586343/...
5920532621959674920092265547231694710826752000,...
-2609410230536249996829038272149121/...
450025282909674287024343687080548396992000,...
307444630821656723644151024596157/...
272742595702832901226874961866999028480,...
-1685658558102333722460472837754189/...
13889669225607231080998261946930506080,...
228262816046878642628376219569959/34211007944845396751227246174705680,...
-2876333034231579924039546769067/18393015024110428360874863534788,1]);
q2 = fliplr([1604612285445093/2658455991569831745807614120560689152,...
2976480701836195/2596148429267413814265248164610048,...
5974590438482987/5070602400912917605986812821504,520154651818787/...
618970019642690137449562112,4332878172069797/9671406556917033397649408,...
3408424760479185/18889465931478580854784,7823297585492293/...
147573952589676412928,1482203634412391/144115188075855872,1]);
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);
p2 = A*(X4*((p2(9)*eye(n))*X4+(p2(5)*eye(n)+p2(6)*X1+p2(7)*X2+p2(8)*X3)*eye(n))+...
(p2(1)*eye(n)+p2(2)*X1+p2(3)*X2+p2(4)*X3));
q2 = X4*((q2(9)*eye(n))*X4+(q2(5)*eye(n)+q2(6)*X1+q2(7)*X2+q2(8)*X3)*eye(n))+...
(q2(1)*eye(n)+q2(2)*X1+q2(3)*X2+q2(4)*X3);
elseif m == 20
X1 = A^2; 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];
p2 = fliplr([504755453995739302254967861091843658659225385703072439/...
19060545079818588288159374768330241280586535602550332513982470267073331200000,...
-8739283052502581091809652173945956169815253539783369/...
114148670977473878836743171447659847170838038103667100934138640957440000,...
7098206776968855944478646605493639889495273479386941/...
73428969634632319719542390989722708706387042054990532764650587750400,...
-113535124075801033808719602712883055821671599003497/...
1648104073175508145683904128728603495708248986725888205771264000,...
159444168407032138165237007441577372335666527110947/...
5304242994128072193005668460275965273543789842336191926620160,...
-250413931005302347677105318749990589129438025771967/...
30554395127465853646346016476244039594146254852167004185600,...
89477294557053179299513803506788857306761134861/...
65142408168740093907440765129293961270139550682600640,...
-12973367660334266409984992764341733133099094183655/...
96446954316495750146294243927538003769401056982850392,...
2886168459322903018906421202468694698344454148801/...
413757847775614543742146048595186631357361891818320,...
-755789275315695418647501273540203715920168810921/...
4774129012795552427793992868405999592584944905596,1]);
q2 = fliplr([2290207460112795/44601490397061246283071436545296723011960832,...
3841593531710827/21778071482940061661655974875633165533184,...
3463624789758205/10633823966279326983230456482242756608,...
8728650941242757/20769187434139310514121985316880384,4208410993374357/...
10141204801825835211973625643008,1606756025554677/...
4951760157141521099596496896,7804153743481773/38685626227668133590597632,...
1845833986288687/18889465931478580854784,1293253388138245/...
36893488147419103232,2408831652010703/288230376151711744,1]);
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);
p2 = A*(X5*((p2(11)*eye(n))*X5+(p2(6)*eye(n)+p2(7)*X1+p2(8)*X2+p2(9)*X3+p2(10)*X4)*eye(n))...
+p2(1)*eye(n)+p2(2)*X1+p2(3)*X2+p2(4)*X3+p2(5)*X4);
q2 = X5*((q2(11)*eye(n))*X5+(q2(6)*eye(n)+q2(7)*X1+q2(8)*X2+q2(9)*X3+q2(10)*X4)*eye(n))...
+(q2(1)*eye(n)+q2(2)*X1+q2(3)*X2+q2(4)*X3+q2(5)*X4);
end
C = q\p;
S = q2\p2;