Satellite := module() description "routines for satellite unknowns determination in linear differential systems"; option package; export Testing, Determination, LinSatTesting, LinearlySatellite; local SatDeterminationN, Matrix2List, BuildTrat, CheckOnlyRational, LinExpr, LinearlySatelliteN; Matrix2List := proc(A::Matrix) # by rows local w, h, r, i, j; h, w := op(1, A); r := []; for i from 1 to h do for j from 1 to w do r := [op(r), A[i,j]]; end do; end do; return r; end proc: BuildTrat := proc(A::Matrix, B::Matrix) local wa, ha, wb, hb, T, _t, R, z; ## testing necessary condition : z' - (trA - trB)z = 0 if LinearFunctionalSystems:-RationalSolution( [diff(z(x),x)-(LinearAlgebra:-Trace(A)-LinearAlgebra:-Trace(B))*z(x)], [z(x)])[1] = 0 then return Matrix([0]); end if; ha,wa := op(1, A); hb,wb := op(1, B); T := Matrix(hb, wa, (i,j)->_t[i,j](x)); #B.T-T.A-map(diff, T, x) R := LinearFunctionalSystems:-RationalSolution( [seq(eq, eq in B.T-T.A-map(diff, T, x))], Matrix2List(T)); T := Matrix(op(1,T), R); return eval(T, [seq(parse(cat("_c[",i,"]=_C",i)),i=1..hb*wa)]); end proc: LinExpr := proc(A, s, v, R) local red, str, m, n, T, l, rr, k, i, res; m := LinearAlgebra:-RowDimension(A); red := OreTools:-Consequences:-ReducedSystem(A,s,R); n := LinearAlgebra:-RowDimension(red[1]); str := red[2]; str := [seq([str[i][1], str[i+1][2]-str[i][2]], i = 1..nops(str)-1), [str[-1][1], n+1-str[-1][2]]]; i := 1; T := Matrix(m); for l in str do rr := Vector[row](m); rr[l[1]] := 1; T[i] := rr; i := i+1; for k from 2 to l[2] do rr := map((e, r)->OreTools:-Apply(OrePoly(0,1), e, r), rr, R) + rr.A; T[i] := rr; i := i+1 end do; end do; T := LinearAlgebra:-LinearSolve(LinearAlgebra:-Transpose(T),Vector(m, i->`if`(i=v,1,0))); res := Matrix(nops(s), max(seq(x[2], x in str))); i := 0; n := 0; for l in str do i := i + 1; m := 0; for k from 1 to l[2] do m := m + 1; n := n + 1; res[i,m] := T[n]; end do; end do; return res; end proc: LinSatTesting := proc(A::Matrix, s::set, v::posint, T) local r, c, ws, S1, S2, R; r,c := op(1, A); if r<>c then error "The matrix should be square"; end if; if (s = {}) or (s intersect {seq(i,i=1..r)} <> s) or (s = {seq(i,i=1..r)}) then error "Incorrect set of selected unknowns indicies"; end if; if (v in s) or (v > r) then error "Incorrect index of testing unknown"; end if; R := OreTools:-SetOreRing(x, 'differential'); S1 := OreTools:-Consequences:-ReducedSystem(A,s,R)[1]; if op(1, S1)[1]=r then if nargs>3 then T := LinExpr(A,s,v,R) end if; return true; end if; ws := s union {v}; S2 := OreTools:-Consequences:-ReducedSystem(A,ws,R)[1]; if op(1, S1)[1]=op(1, S2)[1] then if nargs>3 then T := LinExpr(A,s,v,R) end if; return true; end if; return false; end proc: Testing := proc(A::Matrix, s::set, v::posint) local r, c, ws, S1, S2, R, T; if LinSatTesting(A, s, v) then return true; end if; R := OreTools:-SetOreRing(x, 'differential'); S1 := OreTools:-Consequences:-ReducedSystem(A,s,R)[1]; ws := s union {v}; S2 := OreTools:-Consequences:-ReducedSystem(A,ws,R)[1]; T := BuildTrat(Matrix(op(1,S2),S1), S2); if LinearAlgebra:-Determinant(T)<>0 then return true; end if; return FAIL; end proc: # ws - unselected unknown indicies for testing SatDeterminationN := proc(A::Matrix, s::set, ws) local i, r, c, lws, resT, resU; if nargs < 3 then r,c := op(1,A); lws := {seq(i,i=1..r)} minus s; else lws := ws; end if; resT := {}; resU := {}; for i in lws do if Testing(A,s,i)=true then resT := resT union {i}; else resU := resU union {i}; end if; end do; return [resT, resU, {}]; end proc: CheckOnlyRational := proc(A::Matrix) local Y, r, S, A1, A0, V; r := op(1, A)[1]; A1 := Matrix(r, shape=identity); A0 := -A; if r = 1 then Y := Vector[column]([parse("_y1(x)")]); else Y := Vector[column]([seq(parse(x), x in cat("_y", 1..r, "(x)"))]); end if; S := A1.map(diff, Y, x) + A0.Y; S := [seq(S[i], i = 1..r)]; V := Vector(LinearFunctionalSystems[RationalSolution](S, [seq(Y[i], i=1..r)])); return has(V, parse(cat("_c[", r, "]"))); end proc: Determination := proc(M, s::set) local A1, A0, lws, r, c, ord, i, R, resT, resU, resF, ss, ex, lv; if type(M,Matrix) then return SatDeterminationN(M, args[2..]); end if; if not type(M,list) then error "1st arg should be a matrix or a list"; end if; ord := nops(M); if ord = 1 then return SatDeterminationN(M[1], args[2..]); end if; A1 := M[1]; r, c := op(1, A1); if nargs < 3 then lws := {seq(i,i=1..r)} minus s; else lws := ws; end if; if ord>2 then ord:=ord-1; A1:=Matrix(ord*r,(i,j)->`if`(i=j,1,0)); A1[(ord-1)*r+1..ord*r,(ord-1)*r+1..ord*r]:=M[1]; A0:=Matrix(ord*r, (i,j)->`if`(j=r+i,-1,0)); for i from 0 to ord-1 do A0[(ord-1)*r+1..ord*r,i*r+1..(i+1)*r] := M[ord-i+1]; end do; else A1:=M[1]; A0:=M[2]; end if; R := OreTools:-SetOreRing(x, 'differential'); resT := {}; resU := {}; resF := {}; if not type(Extract,procedure) then error "'Extract' is not a procedure"; end if; for i in lws do ss := s union {i}; ex := Extract(A1, A0, ss, R); if nops(ex)=2 then # only S_d with all selected and testing unknowns ss := {seq(x[2], x in ex[2])}; lv := select(x->x[1]=i,ex[2])[1][2]; ss := ss minus {lv}; if Testing(ex[1], ss, lv) = true then resT := resT union {i}; else resU := resU union {i}; end if; else # S_d and S_a ss := select(x->x[1]=i,ex[4]); if ss <> {} then resT := resT union {i}; next; end if; lv := select(x->x[1]=i,ex[2])[1][2]; if ex[2]={} then # s=0 if CheckOnlyRational(OreTools:-Consequences:-ReducedSystem(ex[1],{lv},R)[1]) then resT := resT union {i}; else resF := resF union {i}; end if; next; end if; if LinearAlgebra:-VectorNorm(ex[3][1..-1, lv]) <> 0 then resT := resT union {i}; next; end if; ss := {seq(x[2], x in ex[2])}; lv := select(x->x[1]=i,ex[2])[1][2]; ss := ss minus {lv}; if LinSatTesting(ex[1], ss, lv) = true then resT := resT union {i}; else resU := resU union {i}; end if; end if; end do; return [resT, resU, resF]; end proc: # ws - unselected unknown indecies for testing LinearlySatelliteN := proc(A::Matrix, s::set, ws) local i, r, c, lws, res; if nargs < 3 then r,c := op(1,A); lws := {seq(i,i=1..r)} minus s; else lws := ws; end if; res := {}; for i in lws do if LinSatTesting(A,s,i) then res := res union {i}; end if; end do; return res; end proc: LinearlySatellite := proc(M, s::set) local A1, A0, lws, r, c, ord, i, R, res, ss, ex, lv; if whattype(M)=Matrix then return LinearlySatelliteN(M, args[2..]); end if; if whattype(M)<>list then error "1st arg should be a matrix or a list"; end if; ord := nops(M); if ord = 1 then return LinearlySatelliteN(M[1], args[2..]); end if; A1 := M[1]; r, c := op(1, A1); if nargs < 3 then lws := {seq(i,i=1..r)} minus s; else lws := ws; end if; if ord>2 then ord:=ord-1; A1:=Matrix(ord*r,(i,j)->`if`(i=j,1,0)); A1[(ord-1)*r+1..ord*r,(ord-1)*r+1..ord*r]:=M[1]; A0:=Matrix(ord*r, (i,j)->`if`(j=r+i,-1,0)); for i from 0 to ord-1 do A0[(ord-1)*r+1..ord*r,i*r+1..(i+1)*r] := M[ord-i+1]; end do; else A1:=M[1]; A0:=M[2]; end if; R := OreTools:-SetOreRing(x, 'differential'); res := {}; if not type(Extract,procedure) then error "'Extract' is not a procedure"; end if; for i in lws do ss := s union {i}; ex := Extract(A1, A0, ss, R); if nops(ex)=2 then # only S_d with all selected unknowns and the testing one ss := {seq(x[2], x in ex[2])}; lv := select(x->x[1]=i,ex[2])[1][2]; ss := ss minus {lv}; if LinSatTesting(ex[1], ss, lv) = true then res := res union {i}; end if; else # S_d and S_a ss := select(x->x[1]=i,ex[4]); if ss <> {} then # testing in S_a res := res union {i}; next; end if; lv := select(x->x[1]=i,ex[2])[1][2]; if {seq(x, x in ex[3][.., lv])} <> {0} then res := res union {i}; next; end if; if nops(ex[2])=1 then # only testing in S_d next; end if; ss := {seq(x[2], x in ex[2])}; lv := select(x->x[1]=i,ex[2])[1][2]; ss := ss minus {lv}; if LinSatTesting(ex[1], ss, lv) = true then res := res union {i}; end if; end if; end do; return res; end proc: end module: