The Matrix Function Toolbox: cosm_pade(A,m,sq) - File Exchange

No BSD License

Highlights from
The Matrix Function Toolbox

arnoldi(A,q1,m) ARNOLDI Arnoldi iteration
ascent_seq(A) ASCENT_SEQ Ascent sequence for square (singular) matrix.
cosm(A) COSM Matrix cosine by double angle algorithm.
cosm_pade(A,m,sq) COSM_PADE Evaluate Pade approximation to the matrix cosine.
cosmsinm(A) COSMSINM Matrix cosine and sine by double angle algorithm.
cosmsinm_pade(A,m) COSMSINM_PADE Evaluate Pade approximations to matrix cosine and sine.
expm_cond(A) EXPM_COND Relative condition number of matrix exponential.
expm_frechet_pade(A,E,k) EXPM_FRECHET_PADE Frechet derivative of matrix exponential via Pade approx.
expm_frechet_quad(A,E,theta,r... EXPM_FRECHET_QUAD Frechet derivative of matrix exponential via quadrature.
fab_arnoldi(A,b,fun,m) FAB_ARNOLDI f(A)*b approximated by Arnoldi method.
funm_condest1(A,fun,fun_frech... FUNM_CONDEST1 Estimate of 1-norm condition number of matrix function.
funm_condest_fro(A,fun,fun_fr... FUNM_CONDEST_FRO Estimate of Frobenius norm condition number of matrix function.
funm_ev(A,fun) FUNM_EV Evaluate general matrix function via eigensystem.
funm_simple(A,fun) FUNM_SIMPLE Simplified Schur-Parlett method for function of a matrix.
logm_cond(A) LOGM_COND Relative condition number of matrix logarithm.
logm_frechet_pade(A,E) LOGM_FRECHET_PADE Frechet derivative of matrix logarithm via Pade approx.
logm_iss(A) LOGM_ISS Matrix logarithm by inverse scaling and squaring method.
logm_pade_pf(A,m) LOGM_PADE_PF Evaluate Pade approximant to matrix log by partial fractions.
mft_test(n) MFT_TEST Test the Matrix Function Toolbox.
mft_tolerance(A) MFT_TOLERANCE Convergence tolerance for matrix iterations.
polar_newton(A) POLAR_NEWTON Polar decomposition by scaled Newton iteration.
polar_svd(A) POLAR_SVD Canonical polar decomposition via singular value decomposition.
polyvalm_ps(c,A,s) POLYVALM_PS Evaluate polynomial at matrix argument by Paterson-Stockmeyer alg.
power_binary(A,m) POWER_BINARY Power of matrix by binary powering (repeated squaring).
quasitriang_struct(R) QUASITRIANG_STRUCT Block structure of upper quasitriangular matrix.
riccati_xaxb(A,B) RICCATI_XAXB Solve Riccati equation XAX = B in positive definite matrices.
rootpm_newton(A,p,c) ROOTPM_NEWTON Coupled Newton iteration for matrix pth root.
rootpm_real(A,p) ROOTPM_REAL Pth root of real matrix via real Schur form.
rootpm_schur_newton(A,p) ROOTPM_SCHUR_NEWTON Matrix pth root by Schur-Newton method.
rootpm_sign(A,p) ROOTPM_SIGN Matrix Pth root via matrix sign function.
signm(A) SIGNM Matrix sign decomposition.
signm_newton(A,scal) SIGNM_NEWTON Matrix sign function by Newton iteration.
sqrtm_db(A,scale) SQRTM_DB Matrix square root by Denman-Beavers iteration.
sqrtm_dbp(A,scale) SQRTM_DBP Matrix square root by product form of Denman-Beavers iteration.
sqrtm_newton(A,scal,maxit) SQRTM_NEWTON Matrix square root by Newton iteration (unstable).
sqrtm_newton_full(A, X0) SQRTM_NEWTON_FULL Matrix square root by full Newton method.
sqrtm_pd(A) SQRTM_PD Square root of positive definite matrix via polar decomposition.
sqrtm_pulay(A,D) SQRTM_PULAY Matrix square root by Pulay iteration.
sqrtm_real(A) SQRTM_REAL Square root of real matrix by real Schur method.
sqrtm_triang_min_norm(T) SQRTM_TRIANG_MIN_NORM Estimated min norm square root of triangular matrix.
sylvsol(A,B,C) SYLVSOL Solve Sylvester equation.
Contents.m
readme.m
View all files

from The Matrix Function Toolbox by Nick Higham
A MATLAB toolbox connected with functions of matrices.

cosm_pade(A,m,sq)

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;

Highlights from The Matrix Function Toolbox

Highlights from
The Matrix Function Toolbox