qd_t =
define 300 function calls x1(t)---x300(t)
Show older comments
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!
Accepted Answer
syms x(t) [1 300]
whos
3 Comments
Tony Cheng
on 26 Aug 2025
syms qd_t(t) [24 1]
qd_t = formula(qd_t)
whos qd_t
Tony Cheng
on 26 Aug 2025
More Answers (2)
You shouldn't use the way described here
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:
49 Comments
Tony Cheng
on 24 Aug 2025
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.
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?
,, 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
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
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:
?
Tony Cheng
on 25 Aug 2025
Tony Cheng
on 26 Aug 2025
Edited: Torsten
on 26 Aug 2025
Tony Cheng
on 26 Aug 2025
Tony Cheng
on 26 Aug 2025
Torsten
on 26 Aug 2025
If you don't want to supply your actual code, make an executable toy example that reproduces your problem.
Tony Cheng
on 27 Aug 2025
Tony Cheng
on 28 Aug 2025
That's what I get with your code:
[mass_matrix_handle,F_f_handle,Jacobian_handle] = symbolic_FD();
numeric_FD(mass_matrix_handle,F_f_handle,Jacobian_handle)
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
on 29 Aug 2025
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
on 30 Aug 2025
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?
Tony Cheng
on 30 Aug 2025
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 ?
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
on 30 Aug 2025
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
on 31 Aug 2025
Tony Cheng
on 31 Aug 2025
Paul
on 31 Aug 2025
If possible, please provide a link to the paper that uses Kane's Method that yields a DAE.
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
on 31 Aug 2025
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.
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
on 31 Aug 2025
I didn't really look at the code, and am not familiar with stabilization methods for this type of problem.
Tony Cheng
on 31 Aug 2025
Tony Cheng
on 31 Aug 2025
Tony Cheng
on 1 Sep 2025
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
on 1 Sep 2025
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.
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
on 2 Sep 2025
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
on 4 Sep 2025
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
on 5 Sep 2025
Tony Cheng
on 9 Sep 2025
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.
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".
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
on 1 Oct 2025
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
on 1 Oct 2025
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
on 1 Oct 2025
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.
idris
on 17 Sep 2025
0 votes
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.
1 Comment
Tony Cheng
on 1 Oct 2025
Categories
Find more on Equation Solving in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
