help links subscribe changes refresh edit

	SandBox
SandBoxZtransform ←	↑	→ SandBox{{7*7}}

SandBox[Polynomial Sequences: Matrix]

Edit detail for SandBox[Polynomial Sequences: Matrix] revision 5 of 11

	1 2 3 4 5 6 7 8 9 10 11
Editor: rrogers Time: 2018/08/16 13:50:58 GMT+0
Note:

changed:
--- Computing Polynomial Sequence Arrays directly from Generating Functions
---
--
--

changed:
--- OGF is similiar except the factorial n! in t^n/n! 
-- OGF/Power Series is similiar except the factorial n! in t^n/n! 

added:
-- Little error checking.  Most Functions expect UniPotent of Nilpotent
-- Arrays.

removed:
---

changed:
-SJ_upper(4)
-- H and Pascals matrix
--
Gen_H(dim) ==
  H : Matrix(FPFI) := new(dim,dim,0)
  for i in 1..(dim-1) repeat H(i+1,i):=i
  H
--
--
--Generate Pascal Matrices with a scalar (x) coefficient
-- exp(x*H)
--
Gen_Pascal(x,dim) ==
	H : Matrix(FPFI):=x*Gen_H(dim)
	Exp_Matrix(x*H)
--
--
--H Left handed Psuedo inverse H_PIl
-- H_PIl*H ~ I
-- inverts H except for the last row.
--
Gen_H_PIl(dim) == 
	H_PIl : Matrix(FPFI) := new(dim,dim,0)
	for i in 1..(dim-1) repeat H_PIl(i,i+1):=1/i
	H_PIl
--
--
-- Calculate Roman's integral equation
-- (exp(x*H)-1)/H
-- Which must be done in proper order
-- Note that in order for the last row
-- to be correct the expansion must be done
-- to dim+1 and then is trimmed to dim.
--
Gen_RI(x,dim) ==
	RI := Gen_Pascal(x,dim+1)
	RI := RI-RI^0
	RI := Gen_H_PIl(dim+1)*RI
	RI := subMatrix(RI,1,dim,1,dim)
	RI
--
--
-- Matrix exponential limited for now
-- Returns exp(Mat)
--
Exp_Matrix(Mat) ==
	n:=nrows(Mat)
	Mat_ret :Matrix(FPFI) := Mat^0
	for i in 1..(n) repeat Mat_ret := Mat_ret +(Mat^i)/factorial(i)
	Mat_ret

-- Log of particular matrix
-- Requires A - I be nilpotent and square; 
-- i.e. A is Unipotent/Lower unittriangular
--
Ln_Uni_Matrix(A) ==
	n := nrows(A)
	m := ncols(A)
	if n=m then
			A1 := A-A^0
			A_ret : Matrix(FPFI) := new(n,n,0)
			if (A1^n = A_ret) then
					for i in 1..n repeat A_ret:=A_ret+(-A1)^i/i
				else
					output("Ln_nil_matrix entered with bad matrix")
		else
			output("Ln_nil_matrix entered with non-square matrix")
			output(A)
	-A_ret
--
-- (A^B) A a matrix, B either integer or matrix
-- A must be Unipotent
--
Exp_Matrix_AB(A,B) ==
	Exp_Matrix(Ln_Uni_Matrix(A)*B)


--
-- Find p in P^-1* A* P = J
-- J is really a shift
-- A must be Stricly Lower Triangular 
--
ST_J(A) ==
	n:=nrows(A)
	m:=ncols(A)
	P : Matrix(FPFI) := new(n,n,0)
	if (n=m) and (A^n = P) then
			P(1,1):=1
			for i in 1..(n-1) repeat P(1..n,i+1):= A*P(1..n,i)
		else
			output("Input to ST_H not strictly triangular")
	P
	
--
-- Same as ST_H except the equation is P^-1*A *P=H
--
ST_H(A) ==
	n:=nrows(A)
	m:=ncols(A)
	P : Matrix(	FPFI) := new(n,n,0)
	if (n=m) and (A^n = P) then
			P(1,1):=1
			for i in 1..(n-1) repeat P(1..n,i+1):= A*P(1..n,i)/(i)
		else
			output("Input to ST_H not strictly triangular")
	P
--
--
-- Start of Examples
--
dim:=5
H:=Gen_H(dim)
GP :=Gen_Pascal(1,dim)
Id := GP^0
--
--
-- Appell Sequences
-- f(t)*exp(x*t)
--
-- Euler Polynmials from WikiPedia
-- 2/(e^t -1) exp(x*t)
--
Euler_P:=2*((GP+1)^-1)
--
-- Bernoulli Polynomials from WikiPedia
-- (t/(e^t -1))* exp(x*t)
--
Bern_Wiki_P := Gen_RI(1,dim)^-1
--
-- Bernoulli of order a from Roman
---- (t/(e^t -1))^a * exp(x*t)
--
Bern_Roman_P := Exp_Matrix_AB(Bern_Wiki_P,-a)
--
-- Hermite Polynomials wiki version
--exp(x*t-(t^2)/2) *exp(x*t)
--
Herm_Wiki_P := Exp_Matrix(-H^2/2)
--
Herm_Roman_P := Exp_Matrix(-v*H^2/2)
-- Hermite/Gould Hopper of order 3
-- exp(x*t+y*t^3)
HGH_P := Exp_Matrix(y*H^3)
--
--
-- Associated Sequences requires Jordan transform
-- exp(x*f(t)) -> SJ_H(f(H))
--
-- Lower Factorial Exponential
-- Roman PP 57
-- exp(x*Ln(t+1)
--
Lower_Fact_Exp := ST_H(Ln_Uni_Matrix(H+1))
-- Just to verify
x:='x;
1
x
x*(x-1)
x*(x-1)*(x-2)
x*(x-1)*(x-2)*(x-3)
--
--  Exponential Polynomials 
-- Roman PP 64
-- Also http://mathworld.wolfram.com/BellPolynomial.html
-- exp(x*(exp(t)-1))
--
Exp_P := ST_H(GP-1)
--
-- 
-
-- Scheffer Sequences require Jordan form
--
-- g(t)*exp(x*f(t)) ->g(H)*ST_H(f(H))
-- 
-- Laguerre Polynomials Roman PP 110
-- (1-t)^(-a-1)*exp(x*(t/(t-1)))
Lag_P  := Exp_Matrix_AB(1-H,(-a-1))*ST_H(H*(H-1)^-1)
-- Checking
vl:= series 'vl
vvl:=(1-vl)^(-a-1)*exp(x*(vl/(vl-1)));
coefficient(vvl,3)
-- Note the "series" give "Power Series" to get the
-- Taylor series one has to multiply by factorial(n)
--


-- Meixner Polynomials
-- Acuctually this is interesting since it is
-- Listed as Scheffer (which is Exponential)
-- but I think Ordinary OGF is a better choice
--
--Meix_P := Exp_Matrix_AB(1-H-H*d,x)*Exp_Matrix_AB(1-H,-x-b)

fricas

--
--

-- Functions to facilitate Transforming Polynomial Generating Functions into 
-- Coefficient arrays.
-- These are coefficient arrays (A) and can be rendered into a polynomial list
-- by A*[1,x,x^2,...]^T
-- EGF Validation can always done for f(x,t)
-- By doing a Taylor series expansion indexed by t
-- OGF/Power Series is similiar except the factorial n! in t^n/n! 
-- has to be applied to the polynomial.
-- Little error checking.  Most Functions expect UniPotent of Nilpotent
-- Arrays.
--
--
--)Clear all
FPFI ==> Fraction(Polynomial(Fraction(Integer)))

Type: Void

fricas

--
--
-- Lower strict Jordan
-- When used on the left of variables x it divides by x
-- with 1/x =0
-- When used on the right of the coefficient matrix it shifts
-- the columns to the left by one.
-- Really a horizontal shift of the coefficient array and 
-- Should be applied A*SJ_lower*[x...]^T
--
SJ_lower(dim) ==
        J : Matrix(Polynomial(Fraction(Integer))) := new(dim,dim,0)
        for i in 1..(dim-1) repeat J(i+1,i):=1
        J

Type: Void

fricas

--
-- Upper strict Jordan
-- When used on the left of variables x it multiplies by x
-- With last entry 0
-- When used on the right of coefficient matrix it shifts
-- the columns to the right by one.
--
SJ_upper(dim) ==
        J : Matrix(FPFI) := new(dim,dim,0)
        for i in 1..(dim-1) repeat J(i,i+1):=1
        J

Type: Void

fricas

--
--
-- H and Pascals matrix
--
Gen_H(dim) ==
  H : Matrix(FPFI) := new(dim,dim,0)
  for i in 1..(dim-1) repeat H(i+1,i):=i
  H

Type: Void

fricas

--
--
--Generate Pascal Matrices with a scalar (x) coefficient
-- exp(x*H)
--
Gen_Pascal(x,dim) ==
        H : Matrix(FPFI):=x*Gen_H(dim)
        Exp_Matrix(x*H)

Type: Void

fricas

--
--
--H Left handed Psuedo inverse H_PIl
-- H_PIl*H ~ I
-- inverts H except for the last row.
--
Gen_H_PIl(dim) == 
        H_PIl : Matrix(FPFI) := new(dim,dim,0)
        for i in 1..(dim-1) repeat H_PIl(i,i+1):=1/i
        H_PIl

Type: Void

fricas

--
--
-- Calculate Roman's integral equation
-- (exp(x*H)-1)/H
-- Which must be done in proper order
-- Note that in order for the last row
-- to be correct the expansion must be done
-- to dim+1 and then is trimmed to dim.
--
Gen_RI(x,dim) ==
        RI := Gen_Pascal(x,dim+1)
        RI := RI-RI^0
        RI := Gen_H_PIl(dim+1)*RI
        RI := subMatrix(RI,1,dim,1,dim)
        RI

Type: Void

fricas

--
--
-- Matrix exponential limited for now
-- Returns exp(Mat)
--
Exp_Matrix(Mat) ==
        n:=nrows(Mat)
        Mat_ret :Matrix(FPFI) := Mat^0
        for i in 1..(n) repeat Mat_ret := Mat_ret +(Mat^i)/factorial(i)
        Mat_ret

Type: Void

fricas

-- Log of particular matrix
-- Requires A - I be nilpotent and square; 
-- i.e. A is Unipotent/Lower unittriangular
--
Ln_Uni_Matrix(A) ==
        n := nrows(A)
        m := ncols(A)
        if n=m then
                        A1 := A-A^0
                        A_ret : Matrix(FPFI) := new(n,n,0)
                        if (A1^n = A_ret) then
                                        for i in 1..n repeat A_ret:=A_ret+(-A1)^i/i
                                else
                                        output("Ln_nil_matrix entered with bad matrix")
                else
                        output("Ln_nil_matrix entered with non-square matrix")
                        output(A)
        -A_ret

Type: Void

fricas

--
-- (A^B) A a matrix, B either integer or matrix
-- A must be Unipotent
--
Exp_Matrix_AB(A,B) ==
        Exp_Matrix(Ln_Uni_Matrix(A)*B)

Type: Void

fricas

--
-- Find p in P^-1* A* P = J
-- J is really a shift
-- A must be Stricly Lower Triangular 
--
ST_J(A) ==
        n:=nrows(A)
        m:=ncols(A)
        P : Matrix(FPFI) := new(n,n,0)
        if (n=m) and (A^n = P) then
                        P(1,1):=1
                        for i in 1..(n-1) repeat P(1..n,i+1):= A*P(1..n,i)
                else
                        output("Input to ST_H not strictly triangular")
        P

Type: Void

fricas

--
-- Same as ST_H except the equation is P^-1*A *P=H
--
ST_H(A) ==
        n:=nrows(A)
        m:=ncols(A)
        P : Matrix(      FPFI) := new(n,n,0)
        if (n=m) and (A^n = P) then
                        P(1,1):=1
                        for i in 1..(n-1) repeat P(1..n,i+1):= A*P(1..n,i)/(i)
                else
                        output("Input to ST_H not strictly triangular")
        P

Type: Void

fricas

--
--
-- Start of Examples
--
dim:=5

$\label{eq1}5$

(1)

Type: PositiveInteger?

fricas

H:=Gen_H(dim)

fricas

Compiling function Gen_H with type PositiveInteger -> Matrix(
      Fraction(Polynomial(Fraction(Integer))))

$\label{eq2}\left[ \begin{array}{ccccc} 0 & 0 & 0 & 0 & 0 \ 1 & 0 & 0 & 0 & 0 \ 0 & 2 & 0 & 0 & 0 \ 0 & 0 & 3 & 0 & 0 \ 0 & 0 & 0 & 4 & 0$

(2)

Type: Matrix(Fraction(Polynomial(Fraction(Integer))))

fricas

GP :=Gen_Pascal(1,dim)

fricas

Compiling function Exp_Matrix with type Matrix(Fraction(Polynomial(
      Fraction(Integer)))) -> Matrix(Fraction(Polynomial(Fraction(
      Integer))))

fricas

Compiling function Gen_Pascal with type (PositiveInteger, 
      PositiveInteger) -> Matrix(Fraction(Polynomial(Fraction(Integer))
      ))

$\label{eq3}\left[ \begin{array}{ccccc} 1 & 0 & 0 & 0 & 0 \ 1 & 1 & 0 & 0 & 0 \ 1 & 2 & 1 & 0 & 0 \ 1 & 3 & 3 & 1 & 0 \ 1 & 4 & 6 & 4 & 1$

(3)

Type: Matrix(Fraction(Polynomial(Fraction(Integer))))

fricas

Id := GP^0

$\label{eq4}\left[ \begin{array}{ccccc} 1 & 0 & 0 & 0 & 0 \ 0 & 1 & 0 & 0 & 0 \ 0 & 0 & 1 & 0 & 0 \ 0 & 0 & 0 & 1 & 0 \ 0 & 0 & 0 & 0 & 1$

(4)

Type: Matrix(Fraction(Polynomial(Fraction(Integer))))

fricas

--
--
-- Appell Sequences
-- f(t)*exp(x*t)
--
-- Euler Polynmials from WikiPedia
-- 2/(e^t -1) exp(x*t)
--
Euler_P:=2*((GP+1)^-1)

$\label{eq5}\left[ \begin{array}{ccccc} 1 & 0 & 0 & 0 & 0 \ -{1 \over 2}& 1 & 0 & 0 & 0 \ 0 & - 1 & 1 & 0 & 0 \ {1 \over 4}& 0 & -{3 \over 2}& 1 & 0 \ 0 & 1 & 0 & - 2 & 1$

(5)

Type: SquareMatrix?(5,Fraction(Integer))

fricas

--
-- Bernoulli Polynomials from WikiPedia
-- (t/(e^t -1))* exp(x*t)
--
Bern_Wiki_P := Gen_RI(1,dim)^-1

fricas

Compiling function Gen_H_PIl with type PositiveInteger -> Matrix(
      Fraction(Polynomial(Fraction(Integer))))

fricas

Compiling function Gen_RI with type (PositiveInteger, 
      PositiveInteger) -> Matrix(Fraction(Polynomial(Fraction(Integer))
      ))

$\label{eq6}\left[ \begin{array}{ccccc} 1 & 0 & 0 & 0 & 0 \ -{1 \over 2}& 1 & 0 & 0 & 0 \ {1 \over 6}& - 1 & 1 & 0 & 0 \ 0 &{1 \over 2}& -{3 \over 2}& 1 & 0 \ -{1 \over{30}}& 0 & 1 & - 2 & 1$

(6)

Type: Matrix(Fraction(Polynomial(Fraction(Integer))))

fricas

--
-- Bernoulli of order a from Roman
---- (t/(e^t -1))^a * exp(x*t)
--
Bern_Roman_P := Exp_Matrix_AB(Bern_Wiki_P,-a)

fricas

Compiling function Ln_Uni_Matrix with type Matrix(Fraction(
      Polynomial(Fraction(Integer)))) -> Matrix(Fraction(Polynomial(
      Fraction(Integer))))

fricas

Compiling function Exp_Matrix_AB with type (Matrix(Fraction(
      Polynomial(Fraction(Integer)))), Polynomial(Integer)) -> Matrix(
      Fraction(Polynomial(Fraction(Integer))))

$\label{eq7}\left[ \begin{array}{ccccc} 1 & 0 & 0 & 0 & 0 \ {{1 \over 2}\ a}& 1 & 0 & 0 & 0 \ {{{1 \over 4}\ {{a}^{2}}}+{{1 \over{12}}\ a}}& a & 1 & 0 & 0 \ {{{1 \over 8}\ {{a}^{3}}}+{{1 \over 8}\ {{a}^{2}}}}&{{{3 \over 4}\ {{a}^{2}}}+{{1 \over 4}\ a}}&{{3 \over 2}\ a}& 1 & 0 \ {{{1 \over{16}}\ {{a}^{4}}}+{{1 \over 8}\ {{a}^{3}}}+{{1 \over{4 8}}\ {{a}^{2}}}-{{1 \over{120}}\ a}}&{{{1 \over 2}\ {{a}^{3}}}+{{1 \over 2}\ {{a}^{2}}}}&{{{3 \over 2}\ {{a}^{2}}}+{{1 \over 2}\ a}}&{2 \ a}& 1$

(7)

Type: Matrix(Fraction(Polynomial(Fraction(Integer))))

fricas

--
-- Hermite Polynomials wiki version
--exp(x*t-(t^2)/2) *exp(x*t)
--
Herm_Wiki_P := Exp_Matrix(-H^2/2)

$\label{eq8}\left[ \begin{array}{ccccc} 1 & 0 & 0 & 0 & 0 \ 0 & 1 & 0 & 0 & 0 \ - 1 & 0 & 1 & 0 & 0 \ 0 & - 3 & 0 & 1 & 0 \ 3 & 0 & - 6 & 0 & 1$

(8)

Type: Matrix(Fraction(Polynomial(Fraction(Integer))))

fricas

--
Herm_Roman_P := Exp_Matrix(-v*H^2/2)

$\label{eq9}\left[ \begin{array}{ccccc} 1 & 0 & 0 & 0 & 0 \ 0 & 1 & 0 & 0 & 0 \ - v & 0 & 1 & 0 & 0 \ 0 & -{3 \ v}& 0 & 1 & 0 \ {3 \ {{v}^{2}}}& 0 & -{6 \ v}& 0 & 1$

(9)

Type: Matrix(Fraction(Polynomial(Fraction(Integer))))

fricas

-- Hermite/Gould Hopper of order 3
-- exp(x*t+y*t^3)
HGH_P := Exp_Matrix(y*H^3)

$\label{eq10}\left[ \begin{array}{ccccc} 1 & 0 & 0 & 0 & 0 \ 0 & 1 & 0 & 0 & 0 \ 0 & 0 & 1 & 0 & 0 \ {6 \ y}& 0 & 0 & 1 & 0 \ 0 &{{24}\ y}& 0 & 0 & 1$

(10)

Type: Matrix(Fraction(Polynomial(Fraction(Integer))))

fricas

--
--
-- Associated Sequences requires Jordan transform
-- exp(x*f(t)) -> SJ_H(f(H))
--
-- Lower Factorial Exponential
-- Roman PP 57
-- exp(x*Ln(t+1)
--
Lower_Fact_Exp := ST_H(Ln_Uni_Matrix(H+1))

fricas

Compiling function Ln_Uni_Matrix with type SquareMatrix(5,Integer)
       -> Matrix(Fraction(Polynomial(Fraction(Integer))))

fricas

Compiling function ST_H with type Matrix(Fraction(Polynomial(
      Fraction(Integer)))) -> Matrix(Fraction(Polynomial(Fraction(
      Integer))))

$\label{eq11}\left[ \begin{array}{ccccc} 1 & 0 & 0 & 0 & 0 \ 0 & 1 & 0 & 0 & 0 \ 0 & - 1 & 1 & 0 & 0 \ 0 & 2 & - 3 & 1 & 0 \ 0 & - 6 &{11}& - 6 & 1$

(11)

Type: Matrix(Fraction(Polynomial(Fraction(Integer))))

fricas

-- Just to verify
x:='x;

Type: Variable(x)

fricas

$\label{eq12}1$

(12)

Type: PositiveInteger?

fricas

$\label{eq13}x$

(13)

Type: Variable(x)

fricas

x*(x-1)

$\label{eq14}{{x}^{2}}- x$

(14)

Type: Polynomial(Integer)

fricas

x*(x-1)*(x-2)

$\label{eq15}{{x}^{3}}-{3 \ {{x}^{2}}}+{2 \ x}$

(15)

Type: Polynomial(Integer)

fricas

x*(x-1)*(x-2)*(x-3)

$\label{eq16}{{x}^{4}}-{6 \ {{x}^{3}}}+{{11}\ {{x}^{2}}}-{6 \ x}$

(16)

Type: Polynomial(Integer)

fricas

--
--  Exponential Polynomials 
-- Roman PP 64
-- Also http://mathworld.wolfram.com/BellPolynomial.html
-- exp(x*(exp(t)-1))
--
Exp_P := ST_H(GP-1)

fricas

Compiling function ST_H with type SquareMatrix(5,Integer) -> Matrix(
      Fraction(Polynomial(Fraction(Integer))))

$\label{eq17}\left[ \begin{array}{ccccc} 1 & 0 & 0 & 0 & 0 \ 0 & 1 & 0 & 0 & 0 \ 0 & 1 & 1 & 0 & 0 \ 0 & 1 & 3 & 1 & 0 \ 0 & 1 & 7 & 6 & 1$

(17)

Type: Matrix(Fraction(Polynomial(Fraction(Integer))))

fricas

--
-- 
-

$\label{eq18}-$

(18)

Type: Variable(-)

fricas

-- Scheffer Sequences require Jordan form
--
-- g(t)*exp(x*f(t)) ->g(H)*ST_H(f(H))
-- 
-- Laguerre Polynomials Roman PP 110
-- (1-t)^(-a-1)*exp(x*(t/(t-1)))
Lag_P  := Exp_Matrix_AB(1-H,(-a-1))*ST_H(H*(H-1)^-1)

fricas

Compiling function Exp_Matrix_AB with type (SquareMatrix(5,Integer)
      , Polynomial(Integer)) -> Matrix(Fraction(Polynomial(Fraction(
      Integer))))

$\label{eq19}\left[ \begin{array}{ccccc} 1 & 0 & 0 & 0 & 0 \ {a + 1}& - 1 & 0 & 0 & 0 \ {{{a}^{2}}+{3 \ a}+ 2}&{-{2 \ a}- 4}& 1 & 0 & 0 \ {{{a}^{3}}+{6 \ {{a}^{2}}}+{{11}\ a}+ 6}&{-{3 \ {{a}^{2}}}-{{1 5}\ a}-{18}}&{{3 \ a}+ 9}& - 1 & 0 \ {{{a}^{4}}+{{10}\ {{a}^{3}}}+{{35}\ {{a}^{2}}}+{{50}\ a}+{24}}&{-{4 \ {{a}^{3}}}-{{36}\ {{a}^{2}}}-{{104}\ a}-{96}}&{{6 \ {{a}^{2}}}+{{42}\ a}+{72}}&{-{4 \ a}-{16}}& 1$

(19)

Type: Matrix(Fraction(Polynomial(Fraction(Integer))))

fricas

-- Checking
vl:= series 'vl

$\label{eq20}vl$

(20)

Type: UnivariatePuiseuxSeries?(Expression(Integer),vl,0)

fricas

vvl:=(1-vl)^(-a-1)*exp(x*(vl/(vl-1)));

Type: UnivariatePuiseuxSeries?(Expression(Integer),vl,0)

fricas

coefficient(vvl,3)

$\label{eq21}{\left( \begin{array}{@{}l} \displaystyle -{{x}^{3}}+{{\left({3 \ a}+ 9 \right)}\ {{x}^{2}}}+{{\left(-{3 \ {{a}^{2}}}-{{15}\ a}-{18}\right)}\ x}+{{a}^{3}}+ \ \ \displaystyle {6 \ {{a}^{2}}}+{{11}\ a}+ 6$

(21)

Type: Expression(Integer)