How to represent the DAE (Differential-Algebraic Equation)¶
To compute consistent initial values of a DAE and to determine the index (and the corresponding condition number) you have to realize the DAE as an example class. Such a class object has the following attributes:
example_name - string delivering a description of the DAE
function_u - True or False depending on whether an additional function u is used or not. If the attribute function_u is absent, this corresponds to False. The function u can be used to fix initial values for some components or/and to realize (nonlinear) relations between components. See an example below.
type_of_DAE - string: ‘proper’, ‘standard’ or ‘linear’
param proper: f and d are given for f(d’(x,t),x,t) param standard: f is given for f(x’,x,t) param linear: f is given as A(t)x’+B(t)x - q(t) where A(t) and B(t) are realized as the functions eval_A and eval_B. svdtol - double, which determines the zero level in rank determination. It has to be specified only in sensitive cases.
The method uses the package algopy for algorithmic differentiation (AD), therefore you have to apply the methods from algopy to code the DAEs, i.e., all used functions like sin, cos, etc. and the vector/matrices like array, zeros, etc. have to be imported from algopy, see also examples below.
An example written in all accepted formulations¶
We give an example in standard formulation
from algopy import zeros, dot, sin
class Kronecker(object):
def __init__(self):
self.type_of_DAE='standard'
self.example_name='Kronecker_standard'
def f(self,y,x,t):
""" linear DAE in Kronecker normal form in standard formulation f(x',x,t) """
n=5
f = zeros(n, dtype=x)
f[0]=y[0]+x[0]
f[1]=y[2]+x[1]
f[2]=y[3]+x[2]
f[3]=y[4]+x[3]
f[4]= x[4]-sin(t)
return f
formulated as proper DAE
from algopy import zeros, dot, sin
class Kronecker(object):
def __init__(self):
self.type_of_DAE='proper'
self.example_name='Kronecker_proper'
def f(self,dp,x,t):
""" linear DAE in Kronecker normal form in proper formulation f(d(x,t)',x,t) """
n=5
f = zeros(n, dtype=x)
f[0]=dp[0]+x[0]
f[1]=dp[1]+x[1]
f[2]=dp[2]+x[2]
f[3]=dp[3]+x[3]
f[4]= x[4]-sin(t)
return f
def d(self,x,t):
""" function d(x,t) """
d = zeros(4, dtype=x)
d[0]=x[0]
d[1]=x[2]
d[2]=x[3]
d[3]=x[4]
return d
and using the coefficient matrices for linear DAEs
from algopy import zeros, dot, sin
class Kronecker(object):
def __init__(self):
self.type_of_DAE='linear'
self.example_name='Kronecker_linear'
def eval_A(self,t):
""" matrix factor of x' """
n=5
A = zeros((n,n),dtype=t)
A[0,0]=1.
A[1,2]=1.
A[2,3]=1.
A[3,4]=1.
return A
def eval_B(self,t):
""" matrix factor of x """
n=5
B = zeros((n,n),dtype=t)
B[0,0]=1.
B[1,1]=1.
B[2,2]=1.
B[3,3]=1.
B[4,4]=1.
return B
def f(self,y,x,t):
""" linear DAE in Kronecker normal form """
n=5
f = zeros(n, dtype=x)
q=zeros(n,dtype=t)
q[4]=sin(t)
f=dot(self.eval_A(t), y) + dot(self.eval_B(t), x) - q
return f
A main program¶
After coding the example we compute consistent initial values as well as the index of the DAE and its condition number at the computed point by calling the module consistentValuesIndex with the following input and output:
param:
example - object of example class
TaylorOrder - integer, order of the highest derivative of Taylor coefficients
alpha - reference vector as ndarray of size dim = dimension of the DAE
x0 - ndarray containing initial guess of length dim*TaylorOrder
t0 - timepoint
- augmentation = False
- if False: using the TaylorOrder
- if True: augmentation of number of Taylor coefficients starting with 1 up to TaylorOrder+1
ftol = 1e-8 - (parameter for python module minimize)
- disp = False - (parameter for python module minimize)
If disp = True also results of the iteration of the minimizer are displayed.
maxiter = 100 - (parameter for python module minimize): maximal number of iterations
- info_level = 0 - integer, higher level implies more information, 0 ... 3
If info_level > 0 then two log-files are generated: example_name_array_date.log and example_name_consi_date.log n array contains Jacobian matrices, projectors, ranks, SVD values etc.n consi contains values related to the computation of consistent initial values and the solution in full accuracy (16 digits)
return:
- solution as ndarray of length dim*(TaylorOrder+1) containing values for x(t0), x’(t0),x”(t0)/2!, ..., x (TaylorOrder)(t0)/(TaylorOrder)!
Caution: Only x(t0), x’(t0),x”(t0)/2!, ..., x (s)(t0)/s! are consistent for s = TaylorOrder+1-index (see [1])
index of the DAE in the computed consistent point
cond - condition number of the index computation
computation_robust - If True -> The relation TaylorOrder+1-index-s==0 is fulfilled
We emphasize that the parameters augmentation, ftol, disp, maxiter, and info_level are optional, but disp = True and info_level > 0 shoud be set in order to obtain information.
Here we give an example
from numpy import array,zeros,pi
from DAEinitialization.consistentValues import consistentValuesIndex
from AD_examples.KroneckerTestDoku.example_Kronecker_proper import Kronecker
ex = Kronecker()
n=5
K = 10
t0=pi/4
alpha=array([1.,1.,1.,1.,1.])
print 'Solution by minimization'
print 'alpha = ',alpha
ftola=1e-14
x0=zeros(K*n)
x0[:n]=alpha
solution, index, cond, robust = consistentValuesIndex(ex,K,alpha,x0,t0,ftol=ftola,disp=True,info_level=2)
The terminal output is:
Solution by minimization
alpha = [ 1. 1. 1. 1. 1.]
*********************************************************
consistentValuesIndex: Computation of consistent solution
*********************************************************
derivativeArray: norm(f_Taylor) 2.17e+00 func: sqrt(yty) 0.00e+00
derivativeArray: norm(f_Taylor) 2.17e+00 func: sqrt(yty) 2.50e+00
derivativeArray: norm(f_Taylor) 9.17e-16 func: sqrt(yty) 2.48e+00
derivativeArray: norm(f_Taylor) 7.49e-16 func: sqrt(yty) 2.43e+00
derivativeArray: norm(f_Taylor) 3.49e-15 func: sqrt(yty) 2.43e+00
derivativeArray: norm(f_Taylor) 3.54e-16 func: sqrt(yty) 2.43e+00
derivativeArray: norm(f_Taylor) 1.11e-16 func: sqrt(yty) 2.43e+00
derivativeArray: norm(f_Taylor) 1.73e-18 Optimization terminated successfully. (Exit mode 0)
Current function value: 2.4319156158
Iterations: 7
Function evaluations: 7
Gradient evaluations: 7
minimize success = True status = 0
*********************************************************
consistentValuesIndex: Consistent solution at t = 0.785398163397
*********************************************************
[[ 1.0000000000000333e+00 7.0710678118654757e-01
-7.0710678118654746e-01 -7.0710678118654757e-01
7.0710678118654746e-01]
[ -1.0000000000000333e+00 -7.0710678118654746e-01
-7.0710678118654757e-01 7.0710678118654746e-01
7.0710678118654757e-01]
[ 5.0000000000001665e-01 -3.5355339059327379e-01
3.5355339059327373e-01 3.5355339059327379e-01
-3.5355339059327373e-01]
[ -1.6666666666667221e-01 1.1785113019775791e-01
1.1785113019775793e-01 -1.1785113019775791e-01
-1.1785113019775793e-01]
[ 4.1666666666668052e-02 2.9462782549439480e-02
-2.9462782549439476e-02 -2.9462782549439483e-02
2.9462782549439476e-02]
[ -8.3333333333336108e-03 -5.8925565098878951e-03
-5.8925565098878960e-03 5.8925565098878951e-03
5.8925565098878968e-03]
[ 1.3888888888889349e-03 -9.8209275164798274e-04
9.8209275164798252e-04 9.8209275164798274e-04
-9.8209275164798252e-04]
[ -1.9841269841270500e-04 5.2323344486093417e-16
1.4029896452114038e-04 -1.4029896452114036e-04
-1.4029896452114038e-04]
[ 2.4801587301588125e-05 0.0000000000000000e+00
-6.5404180607616771e-17 -1.7537370565142548e-05
1.7537370565142544e-05]
[ -2.7557319223986807e-06 -1.1047131989294517e-16
0.0000000000000000e+00 7.2671311786240866e-18
1.9485967294602830e-06]
[ 2.7557319223986809e-07 0.0000000000000000e+00
1.1047131989294517e-17 0.0000000000000000e+00
-7.2671311786240868e-19]]
consistentValuesIndex: index = 4 D = 11 s = 7 (D - index == s ?) = True
condition number of the derivative array considered in the final step of index determination = 4.11e+01
This DAE in Kronecker canonical form has index 4 and is extensively discussed in [1]. In particular, it is explained there why only the values of the first s = TaylorOrder+1-index Taylor coefficients are consistent.
Pendulum example without/with u¶
We illustrate the use of the optional function u that prescribes additional requirements by means of a well-known nonlinear example.
'''
Created on 01.06.2015
@author: estevez/lamour
'''
from algopy import zeros
class Example_pendulum(object):
def __init__(self):
#self.type_of_DAE=None
# proper: f and d are given for f(d'(x,t),x,t)
# standard: f is given for f(x',x,t)
# linear: f is given as A(t)x'+B(t)x - q(t)
# also self.eval_A and self.eval_B are available
self.type_of_DAE='standard'
self.example_name='pendulum_index3'
self.function_u=False # function u is not used
#self.function_u=True # function u is used
def f(self,y,x,t):
""" Pendulum-DAE index 3 """
f = zeros(5, dtype=x)
m=1.
l=1.
g=1.#9.81
f[0]=y[0]-x[2]
f[1]=y[1]-x[3]
f[2]=m*y[2]-x[4]*x[0]
f[3]=m*y[3]-x[4]*x[1]+m*g
f[4]=x[0]**2+x[1]**2-l**2
return f
def u(self,x):
"""
u contains requirements for the initial values
"""
ux=zeros(2, dtype=x)
ux[0]=x[0]**2-x[1]**2-0.5 # components x_1(t_0) and x_2(t_0) lie not on a circle only but also on a hyperbola
ux[1]=x[2]-0.6 # x_1'(t_0) = 0.6
return ux
if __name__=="__main__":
from numpy import array,zeros as npzeros
from DAEinitialization.consistentValues import consistentValuesIndex
ex = Example_pendulum()
n = 5
K = 6
t0= 0.
alpha=array([1,1,0,0,0])
print 'Solution by minimization'
ftola=1e-10
X0=npzeros(K*n)
X0[:n]=alpha
solution,index=consistentValuesIndex(ex,K,alpha,X0,t0,ftol=ftola,disp=True,info_level=2)
The terminal output for function_u = False is:
Solution by minimization
*********************************************************
consistentValuesIndex: Computation of consistent solution
*********************************************************
derivativeArray: norm(f_Taylor) 1.41e+00 func: sqrt(yty) 0.00e+00
derivativeArray: norm(f_Taylor) 1.41e+00 func: sqrt(yty) 3.54e-01
derivativeArray: norm(f_Taylor) 3.06e-01 func: sqrt(yty) 4.44e-01
derivativeArray: norm(f_Taylor) 1.20e-01 func: sqrt(yty) 4.10e-01
derivativeArray: norm(f_Taylor) 2.73e-02 func: sqrt(yty) 4.14e-01
derivativeArray: norm(f_Taylor) 3.59e-03 func: sqrt(yty) 4.14e-01
derivativeArray: norm(f_Taylor) 5.56e-05 func: sqrt(yty) 4.14e-01
derivativeArray: norm(f_Taylor) 1.12e-07 func: sqrt(yty) 4.14e-01
derivativeArray: norm(f_Taylor) 2.39e-10 func: sqrt(yty) 4.14e-01
derivativeArray: norm(f_Taylor) 1.33e-14 Optimization terminated successfully. (Exit mode 0)
Current function value: 0.414213562373
Iterations: 8
Function evaluations: 9
Gradient evaluations: 8
minimize success = True status = 0
*********************************************************
consistentValuesIndex: Consistent solution at t = 0.0
*********************************************************
[[ 7.0710678118602932e-01 7.0710678118706960e-01
-1.2693867358124460e-16 1.2693868019868950e-16
7.0710678118706960e-01]
[ -1.2693866696379970e-16 1.2693866696379970e-16
5.0000000000000555e-01 -4.9999999999926437e-01
3.8081607781919606e-16]
[ 2.5000000000000278e-01 -2.4999999999963218e-01
8.9759222807062377e-17 1.7951838936584310e-16
-7.4999999999889655e-01]
[ 2.9919740659960588e-17 5.9839449887057899e-17
-1.1785113019732217e-01 -2.3570226039534570e-01
1.7951770446029492e-16]
[ -2.9462782549330543e-02 -5.8925565098836424e-02
8.4625685331637612e-17 -5.2892310646807700e-18
1.9308748341390171e-03]
[ 1.6925174124019039e-17 -1.0578647417818407e-18
-4.1393599728831464e-02 2.9439733604323927e-02
4.2214641751921480e-16]
[ -6.8989332881385776e-03 4.9066222673873209e-03
5.3574180688863011e-17 3.0966818387086630e-17
0.0000000000000000e+00]]
consistentValuesIndex: index = 3 D = 7 s = 4 (D - index == s ?) = True
condition number of the derivative array considered in the final step of index determination = 8.93e+00
The terminal output for function_u = True is:
Solution by minimization
*********************************************************
consistentValuesIndex: Computation of consistent solution
*********************************************************
consistentValues: A function u with additional conditions is used.
JacobianDerivativeArray:----- Additional condition u maybe not consistent! Different ranks!
derivativeArray: norm(f_Taylor) 1.62e+00 func: sqrt(yty) 0.00e+00
derivativeArray: norm(f_Taylor) 1.62e+00
JacobianDerivativeArray:----- Additional condition u maybe not consistent! Different ranks!
func: sqrt(yty) 9.36e-01
derivativeArray: norm(f_Taylor) 1.21e+00 func: sqrt(yty) 1.22e+00
derivativeArray: norm(f_Taylor) 2.86e+00 func: sqrt(yty) 1.06e+00
derivativeArray: norm(f_Taylor) 9.01e-01 func: sqrt(yty) 1.27e+00
derivativeArray: norm(f_Taylor) 1.79e+00 func: sqrt(yty) 1.15e+00
derivativeArray: norm(f_Taylor) 6.95e-01 func: sqrt(yty) 1.29e+00
derivativeArray: norm(f_Taylor) 1.39e+00 func: sqrt(yty) 1.31e+00
derivativeArray: norm(f_Taylor) 1.08e-01 func: sqrt(yty) 1.31e+00
derivativeArray: norm(f_Taylor) 1.28e-05 func: sqrt(yty) 1.31e+00
derivativeArray: norm(f_Taylor) 7.93e-15 Optimization terminated successfully. (Exit mode 0)
Current function value: 1.30688530194
Iterations: 8
Function evaluations: 10
Gradient evaluations: 8
minimize success = True status = 0
*********************************************************
consistentValuesIndex: Consistent solution at t = 0.0
*********************************************************
[[ 0.8660254037844386 0.4999999999999999 0.6 -1.0392304845413265
-0.9400000000000003]
[ 0.6 -1.0392304845413265 -0.8140638795573726
-1.4700000000000002 -3.1176914536239804]
[-0.4070319397786863 -0.7350000000000001 -1.6320000000000006
-0.2909845356715714 -2.2049999999999996]
[-0.5440000000000003 -0.0969948452238572 -1.13253028804237
0.9428000000000006 -0.2909845356715712]
[-0.2831325720105925 0.2357000000000002 0.0513400000000003
1.1321723308754725 -0.1570792947919475]
[ 0.0102680000000001 0.2264344661750945 0.5098097629144182
0.3850754705208051 -0.0408234714967926]
[ 0.0849682938190697 0.0641792450868008 0.3435710566229888
-0.0628518584161503 0. ]]
consistentValuesIndex: index = 3 D = 7 s = 4 (D - index == s ?) = True
condition number of the derivative array considered in the final step of index determination = 4.68e+01
Caraxis example without/with svdtol¶
To determine the index of the DAE according to [2], Section 4, several rank considerations are made computing the corresponding singular values. In case the default tolerance does not lead to a robust rank-decision, the parameter svdtol should be set correspondingly in the example.
We illustrate this with a well know example from Test Set for IVP Solvers
'''
Created on 21.10.2017
@author: estevez
'''
from algopy import zeros,sin,sqrt
class Example_car_axis(object):
def __init__(self):
self.type_of_DAE='standard'
self.example_name='car_axis_index3'
#self.svdtol=1e-12
def f(self,y,x,t):
""" Car axis DA index 3 from
https://archimede.dm.uniba.it/~testset/src/problems/caraxis.f
"""
f = zeros(10, dtype=x)
epsilon = 1.0e-2
M = 10.0e+0
L = 1.0e+0
L0 = 0.5e+0
r = 0.1e+0
w = 10.0e+0
g = 1.0e+0
yb = r*sin(w*t)
xb = sqrt(L*L-yb*yb)
K = epsilon*epsilon*M/2e+0
for k in range( 4):
f[k]=y[k]-x[k+4]
xl = x[0]
yl = x[1]
xr = x[2]
yr = x[3]
lam1 = x[8]
lam2 = x[9]
Ll = sqrt(xl**2+yl**2)
Lr = sqrt((xr-xb)**2+(yr-yb)**2)
f[4] = -K*y[4]+(L0-Ll)*xl/Ll +lam1*xb+2e0*lam2*(xl-xr)
f[5] = -K*y[5]+(L0-Ll)*yl/Ll +lam1*yb+2e0*lam2*(yl-yr)-M*epsilon*epsilon*g/2e0
f[6] = -K*y[6]+(L0-Lr)*(xr-xb)/Lr -2e0*lam2*(xl-xr)
f[7] = -K*y[7]+(L0-Lr)*(yr-yb)/Lr -2e0*lam2*(yl-yr)-M*epsilon*epsilon*g/2e0
f[8] = xb*xl+yb*yl
f[9] = (xl-xr)**2+(yl-yr)**2-L*L
return f
if __name__=="__main__":
from numpy import array,zeros as npzeros
from DAEinitialization.consistentValues import consistentValuesIndex
ex = Example_car_axis()
n = 10
K = 4
t0= 0.
L = 1e0
L0 = 0.5e0
xr = L
xl = 0
yr = L0
yl = yr
xra = -L0/L
xla = xra
yra = 0e0
yla = 0e0
lam1 = 0e0
lam2 = 0e0
alpha=array([xl,yl,xr,yr,xla,yla,xra,yra,lam1,lam2])
print 'Solution by minimization'
print 'alpha = ',alpha
ftola=1e-10
X0=npzeros(K*n)
X0[:n]=alpha
solution,index,cond,robust=consistentValuesIndex(ex,K,alpha,X0,t0,ftol=ftola,disp=True,info_level=3)
If this example is run without setting svdtol, then the output reads
Solution by minimization
alpha = [ 0. 0.5 1. 0.5 -0.5 0. -0.5 0. 0. 0. ]
*********************************************************
consistentValuesIndex: Computation of consistent solution
*********************************************************
derivativeArray: norm(f_Taylor) 1.87e+01 func: sqrt(yty) 0.00e+00
derivativeArray: norm(f_Taylor) 1.87e+01 func: sqrt(yty) 1.05e+00
derivativeArray: norm(f_Taylor) 3.15e+00 func: sqrt(yty) 1.05e+00
derivativeArray: norm(f_Taylor) 3.08e-03 func: sqrt(yty) 1.05e+00
derivativeArray: norm(f_Taylor) 5.91e-09 func: sqrt(yty) 1.05e+00
derivativeArray: norm(f_Taylor) 2.11e-15 Optimization terminated successfully. (Exit mode 0)
Current function value: 1.04955261219
Iterations: 4
Function evaluations: 5
Gradient evaluations: 4
minimize success = True status = 0
*********************************************************
consistentValuesIndex: Consistent solution at t = 0.0
*********************************************************
[[ -8.9255898823304907e-29 4.9972575748955489e-01
9.9999999999749289e-01 4.9972351822825956e-01
-4.9972575748955489e-01 2.0140803735178906e-05
-4.9972340731243758e-01 1.0495524674767076e+00
-5.5079932895620248e-04 -2.7538959407692406e-04]
[ -4.9972575748955489e-01 2.0140803735178906e-05
-4.9972340731243758e-01 1.0495524674767076e+00
-4.0281607470357812e-05 -4.5151744578685926e-01
-1.1015583763077088e+00 -4.4703398984210285e-01
5.0371508348089691e-02 1.2591645862606341e-02]
[ -2.0140803735178906e-05 -2.2575872289342963e-01
-5.5077918815385440e-01 -2.2351699492105143e-01
2.4913975406923708e+01 7.1768221052274511e-03
2.4906808720595549e+01 -5.0130584434012647e+01
-4.9757597289641298e-02 -2.4882911261398414e-02]
[ 8.3046584689745693e+00 2.3922740350758170e-03
8.3022695735318504e+00 -1.6710194811337548e+01
2.3993005346636511e-04 -1.9954950512992808e-01
5.9574034153676236e-02 -3.5559068210510442e-02
-3.4563133947197483e-02 -8.6027175253660972e-03]
[ 5.9982513366591278e-05 -4.9887376282482021e-02
1.4893508538419059e-02 -8.8897670526276104e-03
6.9126267894380645e-05 -2.4038571450451165e-04
8.1606318260613912e-05 -3.3433356980115825e-01
0.0000000000000000e+00 0.0000000000000000e+00]]
Check_s_full: rank = 44 used tolerance = 8.96812182213e-13 condition number of B^k = 4.83e+07
*********************************************************
Check_s_full: Determination of s-fullness considering V
*********************************************************
Check_s_full: j = 0 max_V = 4.79935658274e-16 eps*cond = 1.07194398549e-08
Check_s_full: j = 1 max_V = 1.10813798951e-13 eps*cond = 1.07194398549e-08
Check_s_full: j = 2 max_V = 0.0013372490069 eps*cond = 1.07194398549e-08
**************************************************************************************
CheckOmegaTaylorAndProjectors: Projector based decomposition of the derivative array.
**************************************************************************************
A : noOfColumns= 10 Rank= 8 S(1)= 1.000000e+00 S(r)= 1.000000e+00 S(r+1)= 0.000000e+00 S(end)= 0.000000e+00 tol = 6.280370e-15
Computation WR: noOfColumns= 10 Rank= 8 S(1)= 1.000000e+00 S(r)= 5.000000e-04 S(r+1)= 0.000000e+00 S(end)= 0.000000e+00 tol = 4.440893e-15
WR*GL : noOfColumns= 2 Rank= 2 S(1)= 2.920810e+00 S(r)= 6.847416e-01 Full rank. tol = 1.332268e-15
Computation T : noOfColumns= 10 Rank= 8 S(1)= 3.087253e+00 S(r)= 1.000000e+00 S(r+1)= 0.000000e+00 S(end)= 0.000000e+00 tol = 9.155134e-15
Computation WR: noOfColumns= 20 Rank= 16 S(1)= 3.087253e+00 S(r)= 4.472136e-04 S(r+1)= 1.399571e-16 S(end)= 0.000000e+00 tol = 2.945829e-14
WR*GL : noOfColumns= 4 Rank= 4 S(1)= 2.920810e+00 S(r)= 6.847416e-01 Full rank. tol = 3.014092e-15
Computation T : noOfColumns= 10 Rank= 10 S(1)= 3.087253e+00 S(r)= 2.201829e-14 Full rank. tol = 9.809318e-15
Condition CM : noOfColumns= 30 Rank= 24 S(1)= 4.520801e+00 S(r)= 3.329065e-04 S(r+1)= 1.436539e-16 S(end)= 0.000000e+00 tol = 6.649438e-14
consistentValuesIndex: index = 2 D = 5 s = 2 (D - index == s ?) = False
condition number of the derivative array considered in the final step of index determination = 1.36e+04
consistentValuesIndex: The relation D-mu == s is not fulfilled. Check the condition number and the validity of the index determination!!!!!
Check the singular value S(r) in the terminal output or in the array-file and maybe adapt the parameter svdtol in the example.
Note that, apart from this index-determination, we also check the s-fullness according to Corollaries 3 and 4 from [1]. If for the index mu and D=K+1 the relationship D-mu == s is not fulfilled, then some rank determinations may be critical.
A clear hint that the rank-decision was not robust can be found if info_level=3 in the terminal output or in the related array-file (here e.g. car_axis_index3_array_2017-10-26.log)
**************************************************************************************
CheckOmegaTaylorAndProjectors: Projector based decomposition of the derivative array.
**************************************************************************************
A : noOfColumns= 10 Rank= 8 S(1)= 1.000000e+00 S(r)= 1.000000e+00 S(r+1)= 0.000000e+00 S(end)= 0.000000e+00 tol = 6.280370e-15
Computation WR: noOfColumns= 10 Rank= 8 S(1)= 1.000000e+00 S(r)= 5.000000e-04 S(r+1)= 0.000000e+00 S(end)= 0.000000e+00 tol = 4.440893e-15
WR*GL : noOfColumns= 2 Rank= 2 S(1)= 2.920810e+00 S(r)= 6.847416e-01 Full rank. tol = 1.332268e-15
Computation T : noOfColumns= 10 Rank= 8 S(1)= 3.087253e+00 S(r)= 1.000000e+00 S(r+1)= 0.000000e+00 S(end)= 0.000000e+00 tol = 9.155134e-15
Computation WR: noOfColumns= 20 Rank= 16 S(1)= 3.087253e+00 S(r)= 4.472136e-04 S(r+1)= 1.399571e-16 S(end)= 0.000000e+00 tol = 2.945829e-14
WR*GL : noOfColumns= 4 Rank= 4 S(1)= 2.920810e+00 S(r)= 6.847416e-01 Full rank. tol = 3.014092e-15
Computation T : noOfColumns= 10 Rank= 10 S(1)= 3.087253e+00 S(r)= 2.201829e-14 Full rank. tol = 9.809318e-15
Compare the latter singular value, which was considered to be nonzero (S(r) = 2.201829e-14), with the tolerance tol = 9.809318e-15.
In contrast, if svdtol=1e-12, then a robust index characterization is obtained:
**************************************************************************************
CheckOmegaTaylorAndProjectors: Projector based decomposition of the derivative array.
**************************************************************************************
A : noOfColumns= 10 Rank= 8 S(1)= 1.000000e+00 S(r)= 1.000000e+00 S(r+1)= 0.000000e+00 S(end)= 0.000000e+00 tol = 2.828427e-11
Computation WR: noOfColumns= 10 Rank= 8 S(1)= 1.000000e+00 S(r)= 5.000000e-04 S(r+1)= 0.000000e+00 S(end)= 0.000000e+00 tol = 2.000000e-11
WR*GL : noOfColumns= 2 Rank= 2 S(1)= 2.920810e+00 S(r)= 6.847416e-01 Full rank. tol = 6.000000e-12
Computation T : noOfColumns= 10 Rank= 8 S(1)= 3.087253e+00 S(r)= 1.000000e+00 S(r+1)= 0.000000e+00 S(end)= 0.000000e+00 tol = 4.123106e-11
Computation WR: noOfColumns= 20 Rank= 16 S(1)= 3.087253e+00 S(r)= 4.472136e-04 S(r+1)= 1.399571e-16 S(end)= 0.000000e+00 tol = 1.326683e-10
WR*GL : noOfColumns= 4 Rank= 4 S(1)= 2.920810e+00 S(r)= 6.847416e-01 Full rank. tol = 1.357426e-11
Computation T : noOfColumns= 10 Rank= 8 S(1)= 3.087253e+00 S(r)= 1.000000e+00 S(r+1)= 3.852045e-13 S(end)= 2.201829e-14 tol = 4.417724e-11
Computation WR: noOfColumns= 30 Rank= 24 S(1)= 4.668001e+00 S(r)= 1.769493e-06 S(r+1)= 3.159745e-16 S(end)= 3.159745e-16 tol = 3.389377e-10
WR*GL : noOfColumns= 6 Rank= 6 S(1)= 2.921356e+00 S(r)= 6.847196e-01 Full rank. tol = 2.717694e-11
Computation T : noOfColumns= 10 Rank= 10 S(1)= 3.087254e+00 S(r)= 6.847398e-01 Full rank. tol = 5.340063e-11
In this case, the index 3 is correctly determined.