define 300 function calls x1(t)---x300(t)

Hi there,

I am using ode15s to resolve a DAE system. Since there are 300 variables x1(t)~x300(t) in my case, so I feel it is a little bit inconvenient to define them by writing

syms x1(t) … x300(t)

one by one. Are there any efficient codes to define them?

Many thanks!

Answer 1

[Walter Roberson]

syms x(t) [1 300]
whos

  Name      Size            Bytes  Class     Attributes

  t         1x1                 8  sym                 
  x         1x1                 8  symfun              
  x1        1x1                 8  symfun              
  x10       1x1                 8  symfun              
  x100      1x1                 8  symfun              
  x101      1x1                 8  symfun              
  x102      1x1                 8  symfun              
  x103      1x1                 8  symfun              
  x104      1x1                 8  symfun              
  x105      1x1                 8  symfun              
  x106      1x1                 8  symfun              
  x107      1x1                 8  symfun              
  x108      1x1                 8  symfun              
  x109      1x1                 8  symfun              
  x11       1x1                 8  symfun              
  x110      1x1                 8  symfun              
  x111      1x1                 8  symfun              
  x112      1x1                 8  symfun              
  x113      1x1                 8  symfun              
  x114      1x1                 8  symfun              
  x115      1x1                 8  symfun              
  x116      1x1                 8  symfun              
  x117      1x1                 8  symfun              
  x118      1x1                 8  symfun              
  x119      1x1                 8  symfun              
  x12       1x1                 8  symfun              
  x120      1x1                 8  symfun              
  x121      1x1                 8  symfun              
  x122      1x1                 8  symfun              
  x123      1x1                 8  symfun              
  x124      1x1                 8  symfun              
  x125      1x1                 8  symfun              
  x126      1x1                 8  symfun              
  x127      1x1                 8  symfun              
  x128      1x1                 8  symfun              
  x129      1x1                 8  symfun              
  x13       1x1                 8  symfun              
  x130      1x1                 8  symfun              
  x131      1x1                 8  symfun              
  x132      1x1                 8  symfun              
  x133      1x1                 8  symfun              
  x134      1x1                 8  symfun              
  x135      1x1                 8  symfun              
  x136      1x1                 8  symfun              
  x137      1x1                 8  symfun              
  x138      1x1                 8  symfun              
  x139      1x1                 8  symfun              
  x14       1x1                 8  symfun              
  x140      1x1                 8  symfun              
  x141      1x1                 8  symfun              
  x142      1x1                 8  symfun              
  x143      1x1                 8  symfun              
  x144      1x1                 8  symfun              
  x145      1x1                 8  symfun              
  x146      1x1                 8  symfun              
  x147      1x1                 8  symfun              
  x148      1x1                 8  symfun              
  x149      1x1                 8  symfun              
  x15       1x1                 8  symfun              
  x150      1x1                 8  symfun              
  x151      1x1                 8  symfun              
  x152      1x1                 8  symfun              
  x153      1x1                 8  symfun              
  x154      1x1                 8  symfun              
  x155      1x1                 8  symfun              
  x156      1x1                 8  symfun              
  x157      1x1                 8  symfun              
  x158      1x1                 8  symfun              
  x159      1x1                 8  symfun              
  x16       1x1                 8  symfun              
  x160      1x1                 8  symfun              
  x161      1x1                 8  symfun              
  x162      1x1                 8  symfun              
  x163      1x1                 8  symfun              
  x164      1x1                 8  symfun              
  x165      1x1                 8  symfun              
  x166      1x1                 8  symfun              
  x167      1x1                 8  symfun              
  x168      1x1                 8  symfun              
  x169      1x1                 8  symfun              
  x17       1x1                 8  symfun              
  x170      1x1                 8  symfun              
  x171      1x1                 8  symfun              
  x172      1x1                 8  symfun              
  x173      1x1                 8  symfun              
  x174      1x1                 8  symfun              
  x175      1x1                 8  symfun              
  x176      1x1                 8  symfun              
  x177      1x1                 8  symfun              
  x178      1x1                 8  symfun              
  x179      1x1                 8  symfun              
  x18       1x1                 8  symfun              
  x180      1x1                 8  symfun              
  x181      1x1                 8  symfun              
  x182      1x1                 8  symfun              
  x183      1x1                 8  symfun              
  x184      1x1                 8  symfun              
  x185      1x1                 8  symfun              
  x186      1x1                 8  symfun              
  x187      1x1                 8  symfun              
  x188      1x1                 8  symfun              
  x189      1x1                 8  symfun              
  x19       1x1                 8  symfun              
  x190      1x1                 8  symfun              
  x191      1x1                 8  symfun              
  x192      1x1                 8  symfun              
  x193      1x1                 8  symfun              
  x194      1x1                 8  symfun              
  x195      1x1                 8  symfun              
  x196      1x1                 8  symfun              
  x197      1x1                 8  symfun              
  x198      1x1                 8  symfun              
  x199      1x1                 8  symfun              
  x2        1x1                 8  symfun              
  x20       1x1                 8  symfun              
  x200      1x1                 8  symfun              
  x201      1x1                 8  symfun              
  x202      1x1                 8  symfun              
  x203      1x1                 8  symfun              
  x204      1x1                 8  symfun              
  x205      1x1                 8  symfun              
  x206      1x1                 8  symfun              
  x207      1x1                 8  symfun              
  x208      1x1                 8  symfun              
  x209      1x1                 8  symfun              
  x21       1x1                 8  symfun              
  x210      1x1                 8  symfun              
  x211      1x1                 8  symfun              
  x212      1x1                 8  symfun              
  x213      1x1                 8  symfun              
  x214      1x1                 8  symfun              
  x215      1x1                 8  symfun              
  x216      1x1                 8  symfun              
  x217      1x1                 8  symfun              
  x218      1x1                 8  symfun              
  x219      1x1                 8  symfun              
  x22       1x1                 8  symfun              
  x220      1x1                 8  symfun              
  x221      1x1                 8  symfun              
  x222      1x1                 8  symfun              
  x223      1x1                 8  symfun              
  x224      1x1                 8  symfun              
  x225      1x1                 8  symfun              
  x226      1x1                 8  symfun              
  x227      1x1                 8  symfun              
  x228      1x1                 8  symfun              
  x229      1x1                 8  symfun              
  x23       1x1                 8  symfun              
  x230      1x1                 8  symfun              
  x231      1x1                 8  symfun              
  x232      1x1                 8  symfun              
  x233      1x1                 8  symfun              
  x234      1x1                 8  symfun              
  x235      1x1                 8  symfun              
  x236      1x1                 8  symfun              
  x237      1x1                 8  symfun              
  x238      1x1                 8  symfun              
  x239      1x1                 8  symfun              
  x24       1x1                 8  symfun              
  x240      1x1                 8  symfun              
  x241      1x1                 8  symfun              
  x242      1x1                 8  symfun              
  x243      1x1                 8  symfun              
  x244      1x1                 8  symfun              
  x245      1x1                 8  symfun              
  x246      1x1                 8  symfun              
  x247      1x1                 8  symfun              
  x248      1x1                 8  symfun              
  x249      1x1                 8  symfun              
  x25       1x1                 8  symfun              
  x250      1x1                 8  symfun              
  x251      1x1                 8  symfun              
  x252      1x1                 8  symfun              
  x253      1x1                 8  symfun              
  x254      1x1                 8  symfun              
  x255      1x1                 8  symfun              
  x256      1x1                 8  symfun              
  x257      1x1                 8  symfun              
  x258      1x1                 8  symfun              
  x259      1x1                 8  symfun              
  x26       1x1                 8  symfun              
  x260      1x1                 8  symfun              
  x261      1x1                 8  symfun              
  x262      1x1                 8  symfun              
  x263      1x1                 8  symfun              
  x264      1x1                 8  symfun              
  x265      1x1                 8  symfun              
  x266      1x1                 8  symfun              
  x267      1x1                 8  symfun              
  x268      1x1                 8  symfun              
  x269      1x1                 8  symfun              
  x27       1x1                 8  symfun              
  x270      1x1                 8  symfun              
  x271      1x1                 8  symfun              
  x272      1x1                 8  symfun              
  x273      1x1                 8  symfun              
  x274      1x1                 8  symfun              
  x275      1x1                 8  symfun              
  x276      1x1                 8  symfun              
  x277      1x1                 8  symfun              
  x278      1x1                 8  symfun              
  x279      1x1                 8  symfun              
  x28       1x1                 8  symfun              
  x280      1x1                 8  symfun              
  x281      1x1                 8  symfun              
  x282      1x1                 8  symfun              
  x283      1x1                 8  symfun              
  x284      1x1                 8  symfun              
  x285      1x1                 8  symfun              
  x286      1x1                 8  symfun              
  x287      1x1                 8  symfun              
  x288      1x1                 8  symfun              
  x289      1x1                 8  symfun              
  x29       1x1                 8  symfun              
  x290      1x1                 8  symfun              
  x291      1x1                 8  symfun              
  x292      1x1                 8  symfun              
  x293      1x1                 8  symfun              
  x294      1x1                 8  symfun              
  x295      1x1                 8  symfun              
  x296      1x1                 8  symfun              
  x297      1x1                 8  symfun              
  x298      1x1                 8  symfun              
  x299      1x1                 8  symfun              
  x3        1x1                 8  symfun              
  x30       1x1                 8  symfun              
  x300      1x1                 8  symfun              
  x31       1x1                 8  symfun              
  x32       1x1                 8  symfun              
  x33       1x1                 8  symfun              
  x34       1x1                 8  symfun              
  x35       1x1                 8  symfun              
  x36       1x1                 8  symfun              
  x37       1x1                 8  symfun              
  x38       1x1                 8  symfun              
  x39       1x1                 8  symfun              
  x4        1x1                 8  symfun              
  x40       1x1                 8  symfun              
  x41       1x1                 8  symfun              
  x42       1x1                 8  symfun              
  x43       1x1                 8  symfun              
  x44       1x1                 8  symfun              
  x45       1x1                 8  symfun              
  x46       1x1                 8  symfun              
  x47       1x1                 8  symfun              
  x48       1x1                 8  symfun              
  x49       1x1                 8  symfun              
  x5        1x1                 8  symfun              
  x50       1x1                 8  symfun              
  x51       1x1                 8  symfun              
  x52       1x1                 8  symfun              
  x53       1x1                 8  symfun              
  x54       1x1                 8  symfun              
  x55       1x1                 8  symfun              
  x56       1x1                 8  symfun              
  x57       1x1                 8  symfun              
  x58       1x1                 8  symfun              
  x59       1x1                 8  symfun              
  x6        1x1                 8  symfun              
  x60       1x1                 8  symfun              
  x61       1x1                 8  symfun              
  x62       1x1                 8  symfun              
  x63       1x1                 8  symfun              
  x64       1x1                 8  symfun              
  x65       1x1                 8  symfun              
  x66       1x1                 8  symfun              
  x67       1x1                 8  symfun              
  x68       1x1                 8  symfun              
  x69       1x1                 8  symfun              
  x7        1x1                 8  symfun              
  x70       1x1                 8  symfun              
  x71       1x1                 8  symfun              
  x72       1x1                 8  symfun              
  x73       1x1                 8  symfun              
  x74       1x1                 8  symfun              
  x75       1x1                 8  symfun              
  x76       1x1                 8  symfun              
  x77       1x1                 8  symfun              
  x78       1x1                 8  symfun              
  x79       1x1                 8  symfun              
  x8        1x1                 8  symfun              
  x80       1x1                 8  symfun              
  x81       1x1                 8  symfun              
  x82       1x1                 8  symfun              
  x83       1x1                 8  symfun              
  x84       1x1                 8  symfun              
  x85       1x1                 8  symfun              
  x86       1x1                 8  symfun              
  x87       1x1                 8  symfun              
  x88       1x1                 8  symfun              
  x89       1x1                 8  symfun              
  x9        1x1                 8  symfun              
  x90       1x1                 8  symfun              
  x91       1x1                 8  symfun              
  x92       1x1                 8  symfun              
  x93       1x1                 8  symfun              
  x94       1x1                 8  symfun              
  x95       1x1                 8  symfun              
  x96       1x1                 8  symfun              
  x97       1x1                 8  symfun              
  x98       1x1                 8  symfun              
  x99       1x1                 8  symfun              

[Tony Cheng]

Hi Walter,

Many thanks for your support. Also I have a code

qd_t = sym( zeros( 24 , 1 ) ) ;

for i = 1 : length( qd_t )

qd_t( i ) = str2sym( char( ( 'qd_t' )' , num2str( i )' , ('(t)')')' ) ;

end

which can define a column vector that can be applied in the next lines in my file. So are there any more compact or concise methods to define a column vector?

[Walter Roberson]

syms qd_t(t) [24 1]
qd_t = formula(qd_t)

qd_t = 

whos qd_t

  Name       Size            Bytes  Class    Attributes

  qd_t      24x1                 8  sym                

[Tony Cheng]

COOL MAN! U are so amazing!

Answer 2

[Torsten]

You shouldn't use the way described here

https://uk.mathworks.com/help/symbolic/solve-semilinear-dae-system.html

to define a DAE system with 300 unknown functions.

You should directly create the mass matrix M and the vector f of the right-hand side that constitute your system.

Look at the examples on how to use ode15s:

https://uk.mathworks.com/help/matlab/ref/ode15s.html

[Tony Cheng]

Hi Torsten,

Thanks so much for your nice answer and suggestion!

Frankly, I am trying to solve a very typical problem, i.e., the dynamics of a constrained mechanical system. I can format the model as

where q is a the genialized coordinate vector, and are its first and second time derivatives, respectively, and is the Lagrange multiplier vector.

From many publications, I know it is an index-1 DAE system and many researchers use ode15s to numerically resolve this equation. However, I don’t know how to format it into a suitable form that can be used by ode15s, such as the mass matrix, the f vector on the right side of the equation, and the state variables. Frankly, I have the numerical values of q , ,, and at the time t=0. Could you pls give me any hints to rewrite this equation into a form suitable for ode15s?

Best regards

[Torsten]

The fact that there are specialized codes for Mechanical Systems indicates that rewriting the equations and solving them with a standard solver like "ode15s" won't be that easy.

Look at the section "Mechanical Systems" under

https://www.unige.ch/~hairer/software.html

If I were you, I'd take the code given there to solve your problem.

If you want to stick to MATLAB, I'd look in the research articles to see if someone mentions a form of the equations tractable for ODE15S.

You could try to rewrite the system as

q1' = q2

M*q2' + Phi_q^t*lambda = F(tau,q1)-C(q1,q2)*q2-G(q1)

Phi_q*q2' = -Phi_q'*q2

which is a form you could use with ODE15S, but I'm not sure if you would not run into index problems.

Do you know this book:

https://www.springerprofessional.de/mechanical-simulation-with-matlab/19858798

?

[Tony Cheng]

Hi Torsten,

Frankly, I have seen many articles but no description of the form of the equations tractable for ODE15S can be found.

I have downloaded the book

https://www.springerprofessional.de/mechanical-simulation-with-matlab/19858798

However, there are no cases where ode15s are used.

I will try the form you gave me.

Further I will try to format the eqns implicitly and ode15i will be tested.

[Torsten]

As said: I strongly recommend using a code that is especially designed for the solution of mechanical systems. I only know the FORTRAN code by the well-known experts Hairer/Wanner I linked to, but maybe there are other free codes available.

[Tony Cheng]

Hi Torsten,

I have reformatted the eqns of index 1 like this

F_vtr_tilde = F_vtr - C_mtx * qv_t - G_vtr ;

M_mtx_f = [ zeros( 24 ) , M_mtx , zeros( 24 , 18 ) ;

zeros( 18 , 24 ) , PHI_q , zeros( 18 ) ;

eye( 24 ) , zeros( 24 , 42 ) ] ;

F_f = [ F_vtr_tilde - PHI_q' * lambda_t ;

- ( PHI_q_dot + 2 * 5 * PHI_q ) * qv_t - 5 ^ 2 * PHI ;

qv_t ] ;

vars_1 = [ qd_t ; qv_t ; lambda_t ] ;

eqns_1 = M_mtx_f * diff( vars_1 , t ) - F_f ;

isLowIndexDAE( eqns_1 , vars_1 )

ans =

logical

1

I will test this manner to resolve this difficult maths problem.

[Tony Cheng]

Hi Torsten,

In the odeset of my DAE problem using ode15s, I supplied the Jacobian as a function handle in odeset. However, the statistics show that the number of times that the partial derivative matrix is not zero. What is the reason for it? I think with the supplied Jacobian matrix, this number should be zero.

Thanks in advance!

[Torsten]

With "Jacobian", MATLAB means dF_f / d_vars1 . Is it really this 66x66 matrix that you computed and supplied ?

Please include your code, otherwise I can only guess what the problem is.

[Tony Cheng]

Hi Torsten,

The following are my codes

clear ; clc;  close all ; tic  
q_d    = sym( 'q_d'    , [ 24 , 1 ] , 'real' ) ; 
q_v    = sym( 'q_v'    , [ 24 , 1 ] , 'real' ) ; 
lambda = sym( 'lambda' , [ 18 , 1 ] , 'real' ) ; 

PHI_q =  jacobian( PHI   , q_d ) ;      ALPHA = 5 ;  BETA  = 5 ;

Unrecognized function or variable 'PHI'.

PHI_q_dot  = par_mtx_2_rank_vtr( PHI_q , q_d ) * kron( q_v , eye( 24 ) ) ;

M_mtx_f_sym = [ zeros( 24 )      , M_mtx , zeros( 24 , 18 ) ;
                zeros( 18 , 24 ) , PHI_q , zeros( 18      ) ;
                eye( 24 )        , zeros( 24 , 42       ) ] ;   

F_vtr_f_sym = [ F_vtr - C_mtx * qv - G_vtr - PHI_q' * lambda ;
                - ( PHI_q_dot + 2 * ALPHA * PHI_q ) * qv - BETA ^ 2 * PHI ;
               qv ] ;

F_Jacobian_sym = jacobian( F_vtr_f_sym , [ q_d ; q_v ; lambda ] ) ;

syms    q_dv_lad_t(t) [ 66 1 ] ; q_dv_lad_t = formula( q_dv_lad_t)  ;   
qd_t  = q_dv_lad_t(  1 : 24 )  ; 
qv_t  = q_dv_lad_t( 25 : 48 )  ; lad_t      = q_dv_lad_t( 49 : 66 ) ; 

M_mtx_f      = subs( M_mtx_f_sym    ,   q_d , qd_t         ) ;
F_vtr_f      = subs( F_vtr_f_sym    , [ q_d ; q_v ; lambda ] , q_dv_lad_t ) ;
F_Jacobian_f = subs( F_Jacobian_sym , [ q_d ; q_v ; lambda ] , q_dv_lad_t ) ;

eqns_1 = M_mtx_f * diff( q_dv_lad_t , t ) - F_vtr_f  ;   isLowIndexDAE( eqns_1 , q_dv_lad_t ) 

mass_matrix_handle = odeFunction( M_mtx_f      , q_dv_lad_t , 'Sparse'       , true   ) ;   %%% eval( M_mtx_ode ) 
F_f_handle         = odeFunction( F_vtr_f      , q_dv_lad_t , external_force , torque ) ;
Jacobian_handle = odeFunction( F_Jacobian_f , q_dv_lad_t , external_force , torque ) ;   %%% eval( M_mtx_ode ) 

option_dae_1 = odeset( 'Mass'             ,  mass_matrix_handle , ...
    'RelTol'           ,  1 * 10 ^ ( - 6 )  , ...
    'AbsTol'           ,  1 * 10 ^ ( - 7 )  , ...
    'MStateDependence' , 'strong'           , ...
    'MassSingular'     , 'yes'              , ...
    'Stats'            , 'on'               , ...
    'BDF'              , 'off'              , ...
    'Refine'           ,  5                 , ...
    'InitialStep',      t_h * 10^( -5 ) ) ;          

q_all_0 = q_dv_lambda_val(  : , 1 ) ;

for t1 =  0 : t_h : t_h * ( LEN - 1 )  
    
    cnt = round ( t1 / t_h ) + 1 ;  fprintf( 'cnt: %i\n', cnt )  ;  
    
    Jacobian_handle_2 = @( t , q_dv_lad_t ) Jacobian_handle( t , q_dv_lad_t , ( Ext_force( cnt , : ) )' , torque_LA( : , cnt ) ) ;
    F_f_handle_2      = @( t , q_dv_lad_t )      F_f_handle( t , q_dv_lad_t , ( Ext_force( cnt , : ) )' , torque_LA( : , cnt ) ) ;
    
    option_dae_2 = odeset( 'InitialSlope'   ,  qr_va_val_dot_lambda_0( : , cnt ) , ...
        'Jacobian'       , Jacobian_handle_2 ) ;          
    
    option_dae   = odeset( option_dae_1 , option_dae_2 ) ;   
    
    [ tspan , q_val ] = ode15s( F_f_handle_2 , [ cnt - 1 , cnt ] * t_h , q_all_0 , option_dae ) ; 
    
end

Many thanks in advance!

[Torsten]

PHI is undefined, par_mtx_2_rank_vtr is unknown.

[Tony Cheng]

PHI is the 18 * 1 algebraic constraint equations, and par_mtx_2_rank_vtr is a small function I wrote by myself to compute the derivative of a symbolic matrix with respect a rank vector.

Thus, PHI is the function of the symbolic vector q_d,

and PHI_q_dot is the functio of both the symbolic vector q_d and q_v.

[Torsten]

If you don't want to supply your actual code, make an executable toy example that reproduces your problem.

[Tony Cheng]

Hi Torsten, please see the attached file. Pls firstly run I_symbolic_FD.m to get the symbolic results and then I_numeric_FD.m to get the numerical results.

[Tony Cheng]

Dear Torsten, have you run the provided code?

[Torsten]

That's what I get with your code:

[mass_matrix_handle,F_f_handle,Jacobian_handle] = symbolic_FD();

Elapsed time is 40.867402 seconds.

numeric_FD(mass_matrix_handle,F_f_handle,Jacobian_handle)

cnt: 1
23 successful steps
1 failed attempts
170 function evaluations
2 partial derivatives
8 LU decompositions
32 solutions of linear systems
Elapsed time is 2.309397 seconds.
cnt: 2
20 successful steps
0 failed attempts
96 function evaluations
1 partial derivatives
6 LU decompositions
26 solutions of linear systems
Elapsed time is 0.144984 seconds.
cnt: 3
19 successful steps
0 failed attempts
95 function evaluations
1 partial derivatives
6 LU decompositions
25 solutions of linear systems
Elapsed time is 0.115081 seconds.
cnt: 4
18 successful steps
0 failed attempts
94 function evaluations
1 partial derivatives
6 LU decompositions
24 solutions of linear systems
Elapsed time is 0.118874 seconds.
cnt: 5
19 successful steps
1 failed attempts
96 function evaluations
1 partial derivatives
7 LU decompositions
27 solutions of linear systems
Elapsed time is 0.139268 seconds.
cnt: 6
22 successful steps
1 failed attempts
99 function evaluations
1 partial derivatives
7 LU decompositions
30 solutions of linear systems
Elapsed time is 0.103593 seconds.
cnt: 7
20 successful steps
1 failed attempts
100 function evaluations
1 partial derivatives
7 LU decompositions
28 solutions of linear systems
Elapsed time is 0.094056 seconds.
cnt: 8
20 successful steps
0 failed attempts
98 function evaluations
1 partial derivatives
6 LU decompositions
26 solutions of linear systems
Elapsed time is 0.083510 seconds.
cnt: 9
19 successful steps
0 failed attempts
97 function evaluations
1 partial derivatives
6 LU decompositions
25 solutions of linear systems
Elapsed time is 0.083450 seconds.
cnt: 10
23 successful steps
0 failed attempts
101 function evaluations
1 partial derivatives
6 LU decompositions
29 solutions of linear systems
Elapsed time is 0.085572 seconds.
cnt: 11
46 successful steps
12 failed attempts
719 function evaluations
9 partial derivatives
21 LU decompositions
111 solutions of linear systems
Elapsed time is 0.559771 seconds.
cnt: 12
44 successful steps
19 failed attempts
1467 function evaluations
20 partial derivatives
29 LU decompositions
122 solutions of linear systems
Elapsed time is 0.890452 seconds.
cnt: 13
47 successful steps
13 failed attempts
1060 function evaluations
14 partial derivatives
23 LU decompositions
118 solutions of linear systems
Elapsed time is 0.562474 seconds.
cnt: 14
26 successful steps
9 failed attempts
814 function evaluations
11 partial derivatives
16 LU decompositions
72 solutions of linear systems
Elapsed time is 0.359789 seconds.
cnt: 15
44 successful steps
17 failed attempts
1266 function evaluations
17 partial derivatives
26 LU decompositions
122 solutions of linear systems
Elapsed time is 0.438848 seconds.
cnt: 16
30 successful steps
3 failed attempts
405 function evaluations
5 partial derivatives
11 LU decompositions
66 solutions of linear systems
Elapsed time is 0.116305 seconds.
cnt: 17
23 successful steps
1 failed attempts
176 function evaluations
2 partial derivatives
8 LU decompositions
37 solutions of linear systems
Elapsed time is 0.049336 seconds.
cnt: 18
21 successful steps
0 failed attempts
103 function evaluations
1 partial derivatives
6 LU decompositions
31 solutions of linear systems
Elapsed time is 0.030169 seconds.
Elapsed time is 0.030273 seconds.

function [mass_matrix_handle,F_f_handle,Jacobian_handle] = symbolic_FD
%clear ; clc;  close all ; 
tic  
load( 'CONST1_6_in_frame_of_the_EE' ) ;

q1_dv  = sym( 'qr1_dv' , [ 6  , 1 ] , 'real' ) ;  XEE_dv     = sym( 'XEE_dv'     , [ 12 , 1 ] , 'real' ) ;  
q2_dv  = sym( 'qr2_dv' , [ 6  , 1 ] , 'real' ) ;  Bite_force = sym( 'Bite_force' , [  3 , 1 ] , 'real' ) ;  
q3_dv  = sym( 'qr3_dv' , [ 6  , 1 ] , 'real' ) ;  torque     = sym( 'torque'     , [  6 , 1 ] , 'real' ) ;  
q4_dv  = sym( 'qr4_dv' , [ 6  , 1 ] , 'real' ) ;  rho        = 7.8                    ;   %  mass density   :  g/mm^3    7.8 
q5_dv  = sym( 'qr5_dv' , [ 6  , 1 ] , 'real' ) ;  dmter_SM   = 5                      ;   %  diameter of SM : m
q6_dv  = sym( 'qr6_dv' , [ 6  , 1 ] , 'real' ) ;  A_area     = pi  * dmter_SM ^ 2 / 4 ;   %  area of SM     : m^2
lambda = sym( 'lambda' , [ 18 , 1 ] , 'real' ) ; 

[ M_GSM1 , C_GSM1 , G_GSM1 , F_GSM1 , PHI1 ] = A1_MCGF_GSM( q1_dv  , XEE_dv , CONST1 , A_area , rho , torque( 1 ) ) ; 
[ M_GSM2 , C_GSM2 , G_GSM2 , F_GSM2 , PHI2 ] = A1_MCGF_GSM( q2_dv  , XEE_dv , CONST2 , A_area , rho , torque( 2 ) ) ; 
[ M_GSM3 , C_GSM3 , G_GSM3 , F_GSM3 , PHI3 ] = A1_MCGF_GSM( q3_dv  , XEE_dv , CONST3 , A_area , rho , torque( 3 ) ) ; 
[ M_GSM4 , C_GSM4 , G_GSM4 , F_GSM4 , PHI4 ] = A1_MCGF_GSM( q4_dv  , XEE_dv , CONST4 , A_area , rho , torque( 4 ) ) ; 
[ M_GSM5 , C_GSM5 , G_GSM5 , F_GSM5 , PHI5 ] = A1_MCGF_GSM( q5_dv  , XEE_dv , CONST5 , A_area , rho , torque( 5 ) ) ; 
[ M_GSM6 , C_GSM6 , G_GSM6 , F_GSM6 , PHI6 ] = A1_MCGF_GSM( q6_dv  , XEE_dv , CONST6 , A_area , rho , torque( 6 ) ) ; 
[ M_EE   , C_EE   , G_EE   , F_EE          ] = A2_MCGF_EE(  XEE_dv , Bite_force                                   ) ; 

M_mtx = blkdiag( M_GSM1 , M_GSM2 , M_GSM3 , M_GSM4 , M_GSM5 , M_GSM6 , M_EE ) ;
C_mtx = blkdiag( C_GSM1 , C_GSM2 , C_GSM3 , C_GSM4 , C_GSM5 , C_GSM6 , C_EE ) ;
G_vtr =        [ G_GSM1 ; G_GSM2 ; G_GSM3 ; G_GSM4 ; G_GSM5 ; G_GSM6 ; G_EE ] ;
F_vtr =        [ F_GSM1 ; F_GSM2 ; F_GSM3 ; F_GSM4 ; F_GSM5 ; F_GSM6 ; F_EE ] ;

PHI = [ PHI1 ; PHI2 ; PHI3 ; PHI4 ; PHI5 ; PHI6 ] ;  

q_d = [ q1_dv( 1 : 3 ) ; q2_dv( 1 : 3 ) ; q3_dv( 1 : 3 ) ; q4_dv( 1 : 3 ) ; q5_dv( 1 : 3 ) ; q6_dv( 1 : 3 ) ; XEE_dv( 1 :  6 ) ] ; 
q_v = [ q1_dv( 4 : 6 ) ; q2_dv( 4 : 6 ) ; q3_dv( 4 : 6 ) ; q4_dv( 4 : 6 ) ; q5_dv( 4 : 6 ) ; q6_dv( 4 : 6 ) ; XEE_dv( 7 : 12 ) ] ; 

PHI_q      =           jacobian( PHI   , q_d )                           ;      ALPHA = 5 ;  BETA  = 5 ;
PHI_q_dot  = par_mtx_2_rank_vtr( PHI_q , q_d ) * kron( q_v , eye( 24 ) ) ;

M_mtx_f_sym = [ zeros( 24 )      , M_mtx , zeros( 24 , 18 ) ;
                zeros( 18 , 24 ) , PHI_q , zeros( 18      ) ;
                eye( 24 )        , zeros( 24 , 42       ) ] ;   %%% F_vtr_tilde = F_vtr - C_mtx * q_v - G_vtr  ;

% % F_vtr = subs( F_vtr ,   q_d , qd_t ) ;
% % C_mtx = subs( C_mtx , [ q_d , q_v  ] , [ qd_t , qv_t ] ) ;
% % G_vtr = subs( G_vtr ,   q_d , qd_t ) ;
% % PHI       = subs( PHI       ,   q_d , qd_t ) ;
% % PHI_q     = subs( PHI_q     ,   q_d , qd_t ) ;
% % PHI_q_dot = subs( PHI_q_dot , [ q_d , q_v  ] , [ qd_t , qv_t ] ) ;

F_vtr_f_sym = [ F_vtr - C_mtx * q_v - G_vtr - PHI_q' * lambda           ;
             - ( PHI_q_dot + 2 * ALPHA * PHI_q ) * q_v - BETA ^ 2 * PHI ;
                 q_v                                                  ] ;

F_Jacobian_sym = jacobian( F_vtr_f_sym , [ q_d ; q_v ; lambda ] ) ;

syms    q_dv_lad_t(t) [ 66 1 ] ; q_dv_lad_t = formula( q_dv_lad_t)  ;   
qd_t  = q_dv_lad_t(  1 : 24 )  ; 
qv_t  = q_dv_lad_t( 25 : 48 )  ; lad_t      = q_dv_lad_t( 49 : 66 ) ; 

M_mtx_f      = subs( M_mtx_f_sym    ,   q_d , qd_t         ) ;
F_vtr_f      = subs( F_vtr_f_sym    , [ q_d ; q_v ; lambda ] , q_dv_lad_t ) ;
F_Jacobian_f = subs( F_Jacobian_sym , [ q_d ; q_v ; lambda ] , q_dv_lad_t ) ;

eqns_1 = M_mtx_f * diff( q_dv_lad_t , t ) - F_vtr_f  ;   
isLowIndexDAE( eqns_1 , q_dv_lad_t ) ;

mass_matrix_handle = odeFunction( M_mtx_f      , q_dv_lad_t , 'Sparse'   , true   ) ;   %%% eval( M_mtx_ode ) 
F_f_handle         = odeFunction( F_vtr_f      , q_dv_lad_t , Bite_force , torque ) ;
Jacobian_handle    = odeFunction( F_Jacobian_f , q_dv_lad_t , Bite_force , torque ) ;   %%% eval( M_mtx_ode ) 

toc
end

function numeric_FD(mass_matrix_handle,F_f_handle,Jacobian_handle)
%clear ; clc ; close all ; 
tic  
%load( 'symbolic_FD_results'                     ) ;   
load( 'torque_LA_no_load q_all_0 initial slope' ) ;   

load   bite_force.txt 
 
[ rank_NO , LEN ] = size( torque_LA ) ;   chew_force = 0 * 1e+6 * bite_force ;   d_v_vio_norm = zeros( LEN , 2 ) ;  

[ mass_EE , Inertia_EE , molar_0 , TL_in_M , TR_in_M ] = A7_end_effector(  ) ;   t_h = 0.1 ;

option_dae_1 = odeset( 'Mass'             ,  mass_matrix_handle , ...
                       'RelTol'           ,  1 * 10 ^ ( - 6 )  , ...
                       'AbsTol'           ,  1 * 10 ^ ( - 7 )  , ...
                       'MStateDependence' , 'strong'           , ...
                       'MassSingular'     , 'yes'              , ...
                       'Stats'            , 'on'               , ...
                       'BDF'              , 'off'              ) ;  %%%%   'InitialStep',      t_h * 10^( -6 )          
%%%         'Refine'           ,  5                 , ...

q_all_0 = q_dv_lambda_val(  : , 1 ) ;

for t1 =  0 : t_h : t_h * ( LEN - 1 )  
    
cnt = round ( t1 / t_h ) + 1 ;  fprintf( 'cnt: %i\n', cnt )  ;  

tic 

Jacobian_handle_2 = @( t , q_dv_lad_t ) Jacobian_handle( t , q_dv_lad_t , ( chew_force( cnt , : ) )' , torque_LA( : , cnt ) ) ;
F_f_handle_2      = @( t , q_dv_lad_t )      F_f_handle( t , q_dv_lad_t , ( chew_force( cnt , : ) )' , torque_LA( : , cnt ) ) ;
               
%option_dae_2 = odeset( 'Jacobian'     , Jacobian_handle_2 , ...
%                       'InitialSlope' , qr_va_val_dot_lambda_0( : , cnt ) ) ; %%% 
%option_dae_2 = odeset( 'Jacobian'     , Jacobian_handle_2);
option_dae_2 = [];

option_dae   = odeset( option_dae_1 , option_dae_2 ) ;   

[ tspan , q_val ] = ode15s( F_f_handle_2 , [ cnt - 1 , cnt ] * t_h , q_dv_lambda_val(  : , cnt ) , option_dae ) ; 

toc

[ q_all_0 , q_dv_lambda_FD( : , cnt ) ] = deal( ( q_val( end , : ) )' ) ; 

diff_DV_lambda(  : , cnt ) = q_dv_lambda_FD( : , cnt ) - q_dv_lambda_val( : , cnt ) ;

end

toc

tm_ln = 0 : t_h : ( cnt - 1 ) * t_h ;  tm_no = 1 : length( tm_ln ) ;
H0 = figure ; 
subplot( 341 ) ; plot( tm_ln , [ q_dv_lambda_FD( 19 , tm_no ) ; q_dv_lambda_val( 19 , tm_no ) ] ) ;   ylabel('X(mm)' ) ;  legend( 'FD' , 'ID' ) ;
subplot( 345 ) ; plot( tm_ln , [ q_dv_lambda_FD( 20 , tm_no ) ; q_dv_lambda_val( 20 , tm_no ) ] ) ;   ylabel('Y(mm)' ) ;  
subplot( 349 ) ; plot( tm_ln , [ q_dv_lambda_FD( 21 , tm_no ) ; q_dv_lambda_val( 21 , tm_no ) ] ) ;   ylabel('Z(mm)' ) ; 

subplot( 342    ) ; plot( tm_ln , [ q_dv_lambda_FD( 22 , tm_no ) ; q_dv_lambda_val( 22 , tm_no ) ] ) ;   ylabel( '\alpha(rad)' ) ;  
subplot( 346    ) ; plot( tm_ln , [ q_dv_lambda_FD( 23 , tm_no ) ; q_dv_lambda_val( 23 , tm_no ) ] ) ;   ylabel( '\beta(rad)'  ) ;  
subplot( 3,4,10 ) ; plot( tm_ln , [ q_dv_lambda_FD( 24 , tm_no ) ; q_dv_lambda_val( 24 , tm_no ) ] ) ;   ylabel( '\gamma(rad)' ) ; 

subplot( 343 )    ; plot( tm_ln , [ q_dv_lambda_FD( 43 , tm_no ) ; q_dv_lambda_val( 43 , tm_no ) ] ) ;   ylabel('d_X(mm/s)' ) ;   
subplot( 347 )    ; plot( tm_ln , [ q_dv_lambda_FD( 44 , tm_no ) ; q_dv_lambda_val( 44 , tm_no ) ] ) ;   ylabel('d_X(mm/s)' ) ;  
subplot( 3,4,11 ) ; plot( tm_ln , [ q_dv_lambda_FD( 45 , tm_no ) ; q_dv_lambda_val( 45 , tm_no ) ] ) ;   ylabel('d_X(mm/s)' ) ; 

subplot( 344 )    ; plot( tm_ln , [ q_dv_lambda_FD( 46 , tm_no ) ; q_dv_lambda_val( 46 , tm_no ) ] ) ;   ylabel('d_\alpha(rad/s)' ) ;   
subplot( 348 )    ; plot( tm_ln , [ q_dv_lambda_FD( 47 , tm_no ) ; q_dv_lambda_val( 47 , tm_no ) ] ) ;   ylabel('d_\beta( rad/s)' ) ;  
subplot( 3,4,12 ) ; plot( tm_ln , [ q_dv_lambda_FD( 48 , tm_no ) ; q_dv_lambda_val( 48 , tm_no ) ] ) ;   ylabel('d_\gamma(rad/s)' ) ; 
end

[Tony Cheng]

yeah, allthough we specify the Jacobian matrix handle, the partial derivatives are still needed in each count.

[Torsten]

yeah, allthough we specify the Jacobian matrix handle, the partial derivatives are still needed in each count.

No. You can supply them, but if not, ode15s will compute them by finite difference approximations. So you can delete the computation of the Jacobian matrix handle.

Why do you think the partial derivatives are needed ?

How did you arrive at the results you compare with (q_dv_lambda_val) ?

[Tony Cheng]

Hi Torsten,

The reason to specify the Jacobian matrix is, our DAEs are very stiff: some components of y are pretty tiny while some others are extremely huge. From the help center of ode15s, it is critically important to provide the Jacobian.

BTW, could u pls delete these files?

[Torsten]

I think in your case it's important to check whether the constraints are fulfilled in the course of the integration, i.e. whether Φ remains zero.

[Tony Cheng]

The constraints are fulfilled in the integration course. The errors in terms of two-norm sum are 10^(-5) mm. It seems ok for the tolerance.

[Torsten]

And do you have an explanation for the differences between the q_dv_lambda_val and the q_dv_lambda_FD results ? Maybe that you keep torque_LA constant in each subinterval of length t_h ? Did you try whether the integration still works for even smaller values of RelTol and AbsTol ?

[Torsten]

Out of interest, I solved the pendulum problem in cartesian coordinates directly according to the equations you listed above:

Maybe it's worth trying to solve your equations this way because you don't need to put lambda as a solution variable and you don't need to work with a mass matrix.

% Solves the pendulum equation in cartesian coordinates
% mx'' = -2*x*lambda, m*y'' = -m*g - 2*y*lambda, PHI(x,y) = x^2+y^2-L^2
%
% Set parameters
m = 1;
g = 9.81;
L = 10;
% Set initial conditions
x0 = 0;
vx0 = 2;
y0 = 10;
vy0 = 0;
Z0 = [x0;vx0;y0;vy0];
% Set integration interval
tspan = [0 1];
% Call solver
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[T,Z] = ode15s(@(t,z)fun(t,z,m,g),tspan,Z0,options);
% Plot results
% Position
figure(1)
plot(Z(:,1),Z(:,3))

% Constraint error
figure(2)
plot(T,Z(:,1).^2+Z(:,3).^2-L^2)

●

% Lagrange parameter
lambda = zeros(size(T));
for i = 1:numel(T)
  [~,lambda(i)] = fun(T(i),Z(i,:),m,g);
end
figure(3)
plot(T,lambda)

function [res,lambda] = fun(t,z,m,g)
  x = z(1);
  vx = z(2);
  y = z(3);
  vy = z(4);
  dPHIdx = 2*x;
  dPHIdy = 2*y;
  A = [m 0 dPHIdx;0 m dPHIdy;dPHIdx,dPHIdy,0];
  b = [0;-m*g;-2*(vx^2+vy^2)];
  acc = A\b;
  res = zeros(4,1);
  res(1) = z(2);
  res(2) = acc(1);
  res(3) = z(4);
  res(4) = acc(2);
  lambda = acc(3);
end

[Paul]

@Tony Cheng

Instead of solving the DAE, have you considered alternatives, such as Kane's Method, to instead formulate the problem as an ODE that can be integrated with an ode solver?

[Tony Cheng]

Hi Paul, many thx for your advice! I will try to use the Kane's method to build the dynamic model and formulate an ODE system.

[Tony Cheng]

However, in this paper, the Kane's method is used, the final model is still of a DAE form.

[Paul]

@Tony Cheng

If possible, please provide a link to the paper that uses Kane's Method that yields a DAE.

[Torsten]

Writing the equations as

one only needs an ODE integrator to solve for position (q) and velocity (q_dot).

This is accomplished by solving the above linear system of equations for q_doubledot (and lambda, but this is not important in the sequel) using backslash (\) and then setting up the ode system as

dq/dt    = q_dot
dq_dot/dt = q_doubledot 

I doubt there is any method that does not need to integrate for q and q_dot in some way.

[Paul]

@Torsten

I agree that the approach you've described might be viable for the specific problem under discussion.

Assuming I understand the setup of the problem under discussion .... two things to keep in mind.

1. The size of the matrix to be inverted increases as the number of constraints in the system increases. The dimension of q in the problem is not clear to me, though the original question suggests it might be on the order of hundreds (6 coordinates for each of 50 elements in the system?). So we might be starting with M quite large to begin with, and then the size of the matrix to be inverted grows with number of constraints (i.e., the dimension of lamda). With Kane's Method the size of the matrix to be inverted decreases as the number of constraints increases, which results in computational savings (which might or might not be significant).

2. The suggested approach will inherently be influenced by computational innaccuracy, which can result in the integrated solution drifting and not satisfying the constraint equation Phi*qdot = 0 (or its integral if the constraints apply at the configuration level). I know that you're aware of this issue as shown in your "constraint error" plot, but I thought it worth mentioning again. With Kane's method, the constraints are inherently incoporated into the EOM for the minimal* number of coordinates and therefore "constraint drift" is not a concern.

*I believe that there are formulations of Kane's Method for a non-minimal number of coordinates, but I'm not familiar with that approach nor its motivations.

@Tony Cheng

Maybe Simscape Multibody would be of interest.

[Torsten]

@Paul

I cannot recommend a certain method because I have no experience with the constrained mechanical system equations. That's why I gave a link to a FORTRAN package that solves this special type of DAE system. Since the authors of the package are experts in the solution of differential-algebraic problems, I guess they took care of the drift problem you described above.

But as @Tony Cheng already prepared everything in MATLAB, I feel it's worth testing the method I used in the solution of the pendulum equations for his problem because I guess that only little changes to his code will be necessary. Further, comparison between the results from his formulation and the new one could be of interest.

Looking again at @Tony Cheng 's code, I'm surprised about his right-hand side of the DAE-system. In particular, I cannot make sense of the line

PHI_q * q_v_dot =  - ( PHI_q_dot + 2 * ALPHA * PHI_q ) * q_v - BETA ^ 2 * PHI ;

Is this some stabilization of the method used ?

[Paul]

I didn't really look at the code, and am not familiar with stabilization methods for this type of problem.

[Tony Cheng]

Hi Torsten, the code incorporates ALPHA and BETA

PHI_q * q_v_dot = - ( PHI_q_dot + 2 * ALPHA * PHI_q ) * q_v - BETA ^ 2 * PHI , because I use the Baumgarte stablization method (well-known in constrained mechanical systems) to try the best to reduce the position, velocity, and acceleration drift errors.

[Tony Cheng]

Hi Paul, the attached paper is the one using Kane's method to generate the EOMs in the form of DAEs.

As far as Simscape Mulitbody is concerned, I will try to do the simuation work in this environment after the coding is finished.

[Tony Cheng]

Hi Torsten, the final form of the EOMs to be implemented in Matlab using ode15s is

the variable vector is [ q; D_q ; lambda ].

[Paul]

@Tony Cheng

Thanks for posting that paper.

I gather that the DAE is shown in equation (31), though the upper-right term on the LHS should be transpose(phi_q).

From what I can tell after a quick skim, the authors intentionally formulate the problem to yield a DAE to avoid "sophisticated velocities transformation between the joint-space and task-space." But that was a choice, not a requirement of the method.

Though the authors are solving a DAE, they explicitly say that they use Matlab ode45. Perhaps they are using the method suggesed by @Torsten.

The coefficient matrix on the LHS in this comment is singular. Is that expected?

[Torsten]

@Tony Cheng

Looking again in your code, there is another thing I don't understand.

In this line

[ tspan , q_val ] = ode15s( F_f_handle_2 , [ cnt - 1 , cnt ] * t_h , q_dv_lambda_val(  : , cnt ) , option_dae ) ;

you choose q_dv_lambda_val as initial condition. Shouldn't this be q_dv_lambda_FD because you want to continue an integration with solution variable q_dv_lambda_FD ? Otherwise, the integrations in the loop are all independent from another since q_dv_lambda_val is prescribed externally - thus there is no temporal continuation of a process.

[Tony Cheng]

At the instant cnt*delta_h, the actuating torques fed into the generalized force vector are updated.

[Torsten]

At the instant cnt*delta_h, the actuating torques fed into the generalized force vector are updated.

That's right, but you must continue the computation with the solution so far obtained with ode15s, not with a solution obtained by an external software (or whereever the q_dv_lambda_val array stems from).

[Tony Cheng]

Yes, u are rite. Now I see the meaning of long time simulations. if q_dv_lambda_val is used as initial condition for each loop, then the integrations in one loop are all independent from another, having nothing to do with the long time simulations.

However, if the q_dv_lambda_FD obtained from the end of the last loop is used as the initial condition of the next loop, the violation in terms of displacement, velocity, and acceleration is very huge in the entire time span.

[Torsten]

However, if the q_dv_lambda_FD obtained from the end of the last loop is used as the initial condition of the next loop, the violation in terms of displacement, velocity, and acceleration is very huge in the entire time span.

I noticed this, too. So either the models used to compute q_dv_lambda_FD and q_dv_lambda_val are not comparable or you or the one responsible for the q_dv_lambda_val results made a mistake in the implementation.

[Tony Cheng]

That could be. From the articles I find that the index-1 DAEs have no acceleration violations. Violations only exist in coordinates and velocity levels. However, from my codes, the acceleration violation is the largest among the three.

Probably I will turn to Simscape multibody for some simulations @Paul.

[Torsten]

However, from my codes, the acceleration violation is the largest among the three.

In this case, I don't think the solution method is the problem. I suspect there are discrepancies in M, PHI and/or F between the two solutions.

[Tony Cheng]

These days I am reading relevant papers and find that in the form

the largest contraint violations occur in displacement and velocity levels, while the violations in acceleration are the smallest. I agree with what you suspect.

[Torsten]

Displacement and velocity are the solution variables, acceleration is a deduced quantity that is computed from these solution variables. Thus violations in acceleration only follow from violations in displacement and velocity (and force terms) - they should not occur per se.

This can easily be seen if you look at the underlying system:

Acceleration and Lagrage parameters can directly be computed from displacement and velocity (and force terms) by solving a linear system of equations. Thus they are only "wrong" if q, qdot, F, C and/or G are "wrong".

[Torsten]

Using the above linear system to compute acceleration and Lagrange parameters from displacement, velocity and forces, I suggest computing acceleration and Lagrange parameters with your code using the displacements and velocities saved in "q_dv_lambda_val". They should be similar to the accelerations and Lagrange parameters of the other code if you expect similar results for the two simulations.

[Tony Cheng]

Hi Stephen,

These days I am always thinking the constraint violation issue. I guess the reason for it could be: since the ode15s resolves an initial value problem, I input the generalized force vector F(tau, q) at t0 and it is employed in every step from t0 to tf. However, tau is time-variant and needs to be updated at every time step in [t0, tf] with the numerically computed q and q_dot. For the case in Solve Differential Algebraic Equations (DAEs) the example from Help Center in matlab, it is OK to use ode15s cause the generalized force is the gravity of the pendulum which is constant. While in my problem, it is probably not right to update tau at each time step between t0 and tf.

Here I want to know, when ode15s is applied, are there any manners by which I can update tau between t0 and tf? If not, perhaps I need to write codes from the algorithms of ode15s.

Best regards

Tony

[Torsten]

But don't you update the forces here

Jacobian_handle_2 = @( t , q_dv_lad_t ) Jacobian_handle( t , q_dv_lad_t , ( chew_force( cnt , : ) )' , torque_LA( : , cnt ) ) ;
F_f_handle_2      = @( t , q_dv_lad_t )      F_f_handle( t , q_dv_lad_t , ( chew_force( cnt , : ) )' , torque_LA( : , cnt ) ) ;

?

My interpretation was that you prescribe new values for chew_force and torque at the beginning of each new timestep. What else do you intend to do now ?

By the way: Who is Stephen ?

[Tony Cheng]

hi Torsten, i mean, at a time span from t0 to tf, the tau is updated at t0 only. while at the time steps in the span, tau is not updated.

[Torsten]

But tau is the time, isn't it ? And the time is permanently updated by ode15s.

F(tau,q) is computed from q at time tau and the values for chew force and torque_LA that are kept constant over a time interval [ cnt - 1 , cnt ] * t_h and are inserted in F_vtr and G_vtr.

[Tony Cheng]

tau is not meaning time. it means torque. so torque is prescribed at t0 only, and at every time step in the time span from t0 to tf, torque is kept constant as at at t0. However, this is not right i think, since torque is the function of q and q_dot, which are different at every time step in the time span.

[Torsten]

Ok, I didn't know that tau is torque.

At the moment, you prescribe torque (I guess, differently) in every time interval [ cnt - 1 , cnt ] * t_h from the external simulation and keep it constant in each interval.

If you think that it has to be computed internally from q and q_dot, you have to compute

G_vtr =        [ G_GSM1 ; G_GSM2 ; G_GSM3 ; G_GSM4 ; G_GSM5 ; G_GSM6 ; G_EE ] ;
F_vtr =        [ F_GSM1 ; F_GSM2 ; F_GSM3 ; F_GSM4 ; F_GSM5 ; F_GSM6 ; F_EE ] ;

differently - thus express torque(1),...,torque(6) as functions of q1_dv,...,q6_dv.

But I'm confused: In your matrix equation, the external forces are explicitly written as F(tau,q). If tau is torque, for me this would mean that torque is independently given and not a function of q. But I have no background knowledge in the field of constrained mechanical systems.

Answer 3

[idris]

You do not need to type out all 300 variables manually, for MATLAB supports creating symbolic functions in bulk. It is recommended to use vector of symbolic functions, which creates x(t) as a 300×1 vector:

syms x(t) [300 1]

Now, one can access them as x(1), x(2), …, x(300), which correspond to x1(t),x2(t),…,x300(t)x_1(t), x_2(t), \ldots, x_{300}(t)x1​(t),x2​(t),…,x300​(t); thereby, providing the cleanest approach which works very well when formulating DAEs for ode15s.

[Tony Cheng]

all rite man, many thx for your nice help!