How to fix error in my DAE program?

The follwing is my main code- I am getting error -

Not enough input arguments.

Error in @(t,y,ydot)reduced21(t,y,ydot,T_a)

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode15s (line 152)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

function [yp] = reduced21(t,y,ydot,T_a)

% Define the initial state vector y0 and the initial derivative vector ydot0
y0 = zeros(12,1);
ydot0 = zeros(12,1);
beta = 13*(pi/180);                      %Helix angle (rad)
speed = 1000/60; 

Freq = 1000*20/60;

% Common parameters
theta_a_vec = -speed*2*pi*t;     %torsional angle of driver gear (rad)
% theta_p = 22.9;               %torsional angle of driven gear (rad)
theta_p = pi/6 + speed*2*pi*t;
%e = 20e-6;                      %circumferential relative displacement of the teeth (m)
e = 0;
z = e*tan(beta);                %axial tooth displacement caused by internal excitation (m)
e_a = 48e-6;                    %circumferential displacement of the driver gear (m)
e_p = 48e-6;                    %circumferential displacement of the driver gear (m)
ks = 0.9e8 + 20*sin(2*pi*Freq*t);                   %Shear stiffness 
Cs = 0.1;                                           %Shear damping


% Driver gear  
m_a = 0.5;                      %mass of driver gear (kg)
c_ay=500;                       %Damping of driver gear in y direction (Ns/m)
c_az=500;                       %Damping of driver gear in z direction (Ns/m)
k_ay= 8e7;                      %Stiffness of driver gear in y direction (N/m)
k_az= 5e7;                      %Stiffness of driver gear in z direction (N/m)
R_a = 51.19e-3;                 %Radius
I_a = 0.5*m_a*(R_a^2);          %Moment of inertia of driver gear (Kg.m3)
y_ac= e_a + theta_a_vec*R_a;    %circumferential displacement at the meshing point on the driver gear (m)
y_pc = e_p - theta_a_vec*R_a;   %circumferential displacement at the meshing point on the driven gear (m)
z_a = e_a*tan(beta);  
z_p = e_p*tan(beta);  
z_ac = (z_a-y_ac)*tan(beta);    %axial displacements of the meshing point on the driving gear  (m)
z_pc = (z_p+y_pc)*tan(beta);    %axial displacements of the meshing point on the driven gear   (m)

% Excitation forces 
Fy=ks*cos(beta)*(y_ac-y_pc-e) + Cs*cos(beta)*2*R_a*speed*2*pi;                %axial dynamic excitation force at the meshing point (N)
Fz=ks*sin(beta)*(z_ac-z_pc-z) - Cs*sin(beta)*2*tan(beta)*R_a*speed*2*pi;      %circumferential dynamic excitation force at the meshing point (N)


%Time interpolation - Needed for solution of equations (It basically uses
%your torque and prescribed time matrices to generate a time matrix to be
%used in the ODE solver)
time = 0:0.00001:0.6;
Torque = interp1(time,T_a,t);
Opp_Torque = 68.853;                        % Average torque value
%Driver gear equations
yp = zeros(12,1);                           %vector of 12 equations
yp(1) = y(2); 
yp(2) = (-Fy - k_ay*y(1) - c_ay*y(2))/m_a;  %acceleration in y (up and down) direction (m/s2)
yp(3) = y(4);  
yp(4) = (-Fz - k_az*y(3) - c_az*y(4))/m_a;  %acceleration in z (side to side) direction (m/s2)
yp(5) = y(6);  
yp(6) = (-Torque - Fy*R_a + Opp_Torque)/I_a;                % angular acceleration (rad/s2)

% Driven gear  
m_p = 0.5;                                  %mass of driven gear (kg)
c_py=500;                                   %Damping of driven gear in y direction (Ns/m)
c_pz=500;                                   %Damping of driven gear in z direction (Ns/m)
k_py= 9.5e7;                                %Stiffness of driven gear in y direction (N/m)
k_pz =9e7;                                  %Stiffness of driven gear in z direction (N/m)
R_p = 139.763e-3;                           %Radius
I_p = 0.5*m_p*(R_p^2);                      %Moment of inertia of driven gear (Kg.m3)
n_a = 20;                           % no of driver gear teeth  
n_p = 60;                           % no of driven gear teeth  
i = n_p/n_a;                        % gear ratio

%Driven gear equations
yp(7) = y(8);  
yp(8) = (-Fy - k_py*y(7) - c_py*y(8))/m_p;            %acceleration in y (up and down) direction (m/s2)
yp(9) = y(10);  
yp(10) = (-Fz - k_pz*y(9) - c_pz*y(10))/m_p;          %acceleration in y (side to side) direction (m/s2)
yp(11) = y(12);  
yp(12) = (-Torque*i - Fy*R_p + Opp_Torque*i)/I_p;                       % angular accelration (rad/s2)

y0 = [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0]; % matrix of intial values of y
ydot0 = [0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0; 0]; %matrix of initial vlaues of ydot0

end  

My command window code-

tic
A = load("torque_for_Sid_60s.mat");
T_a = A.torque_60s(1:60001);                  %Torque (Nm)

speed = 1000/60; 
tspan=0:0.00001:0.6;  
%tspan = [0 12];  
y0 = zeros(12,1);
ydot0 = zeros(12,1);
y0 = zeros(12,1);  
[t,y,ydot] = ode15s(@(t,y,ydot) reduced21(t,y,ydot,T_a),tspan,y0,ydot0);
toc

% Driver gear graphs
figure
subplot(3,1,1)
plot(t,y(:,2))
ylabel('(y) up and down acceleration (m/s2)')
xlabel('Time')
title('acceleration in y direction vs time')
grid on

subplot(3,1,2)
plot(t,y(:,4))
ylabel('(z) side to side acceleration (m/s2)')
xlabel('Time')
title('acceleration in z direction vs time')
grid on

subplot(3,1,3)
plot(t,y(:,6))
ylabel('angular acceleration (rad/s2)')
xlabel('Time')
title('angular acceleration vs time')
grid on

% Driven gear graphs
figure
subplot(3,1,1)
plot(t,y(:,8))
ylabel('(y) up and down acceleration (m/s2)')
xlabel('Time')
title('acceleration in y direction vs time')
grid on

subplot(3,1,2)
plot(t,y(:,10))
ylab
el('(z) side to side acceleration (m/s2)')
xlabel('Time')
title('acceleration in z direction vs time')
grid on
 
subplot(3,1,3)
plot(t,y(:,12))
ylabel('angular acceleration (rad/s2)')
xlabel('Time')
title('angular acceleration vs time')
grid on
 
% Torque graph
figure
plot(t,T_a(1:60001))
ylabel('Torque (Nm)')
xlabel('Time (sec)')
title('Torque vs time')
grid on

Answer 1

[Torsten]

The call to ode15s is incorrect:

[t,y] = ode15s(@(t,y) reduced21(t,y,T_a),tspan,y0);

instead of

[t,y,ydot] = ode15s(@(t,y,ydot) reduced21(t,y,ydot,T_a),tspan,y0,ydot0);

and the function header has to be modified accordingly:

function yp =  reduced21(t,y,T_a)
...
end

[Siddharth Jain]

Thanks, but the values of y0 and ydot0 are not getting implemented into the function

[Torsten]

You set

y0 = zeros(12,1);

before you call ode15s. So it will be used as vector of initial values for the unknowns.

yp0 does not need to be supplied because it can be calculated using your function reduced21.

Or do you use ydot somewhere in the function ? I can't see it.

[Siddharth Jain]

ydot0 is a vector of differential of the initial values in y0

[Torsten]

It's not needed (and you cannot prescribe it).

ydot0 will be the vector yp that is calculated in reduced21 by inserting (t,y) = (0,y0) into the equations.