How to convert ode45 to Differential Algebraic Equations?

I have 12 equations with 6 DOF ( see yp) that i want to convert to Differential Algebraic Equation format. How can I do this?

function [yp] = reduced(t,y,T_a)  
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)
e = 20e-6;                      %circumferential relative displacement of the teeth (m)
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 = 5.3*e-2;                  %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: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 = 2.5*e-1;                              %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= 8e7;                                  %Stiffness of driven gear in y direction (N/m)
k_pz =5e7;                                  %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
Opp_Torque = 68.853;                % Average torque value

%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)

end  

[Torsten]

What do you mean by " convert to Differential Algebraic Equation format " ?

[Siddharth Jain]

I want to use DAE for simulataneously solving the equations instead of using ode45

[Torsten]

Solving differential equations and algebaic equations simultaneously is not possible with ODE45. If this is what your question aims at, you will have to use ODE15S or ODE23T, as John D'Errico colourfully explained.

Answer 1

[John D'Errico]

I would like to fly to Los Angeles. However, I don't want to take a plane, and since I have often taken a train to destinations, does anyone know how to make a train fly?

The point is, a train is not designed to fly. No matter how hard you try, you will not get it off the tracks without rather unhappy consequences.

In your case, you have a DAE: a Differential Algebraic Equation system. So you use a tool designed to solve a DAE. As I recall, ODE15S will do that, and ODE23T. Maybe some others, but those are the ones that spring to mind.

help ode15s

 ODE15S Solve stiff differential equations and DAEs, variable order method.
    [TOUT,YOUT] = ODE15S(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates 
    the system of differential equations y' = f(t,y) from time T0 to TFINAL 
    with initial conditions Y0. ODEFUN is a function handle. For a scalar T 
    and a vector Y, ODEFUN(T,Y) must return a column vector corresponding 
    to f(t,y). Each row in the solution array YOUT corresponds to a time
    returned in the column vector TOUT.  To obtain solutions at specific
    times T0,T1,...,TFINAL (all increasing or all decreasing), use TSPAN =
    [T0 T1 ... TFINAL].     
    
    [TOUT,YOUT] = ODE15S(ODEFUN,TSPAN,Y0,OPTIONS) solves as above with default
    integration properties replaced by values in OPTIONS, an argument created
    with the ODESET function. See ODESET for details. Commonly used options
    are scalar relative error tolerance 'RelTol' (1e-3 by default) and vector
    of absolute error tolerances 'AbsTol' (all components 1e-6 by default).  
    If certain components of the solution must be non-negative, use
    ODESET to set the 'NonNegative' property to the indices of these
    components.  The 'NonNegative' property is ignored for problems
    where there is a  mass matrix. 
    
    The Jacobian matrix df/dy is critical to reliability and efficiency. Use
    ODESET to set 'Jacobian' to a function handle FJAC if FJAC(T,Y) returns 
    the Jacobian df/dy or to the matrix df/dy if the Jacobian is constant. 
    If the 'Jacobian' option is not set (the default), df/dy is approximated 
    by finite differences. Set 'Vectorized' 'on' if the ODE function is coded 
    so that ODEFUN(T,[Y1 Y2 ...]) returns [ODEFUN(T,Y1) ODEFUN(T,Y2) ...]. 
    If df/dy is a sparse matrix, set 'JPattern' to the sparsity pattern of
    df/dy, i.e., a sparse matrix S with S(i,j) = 1 if component i of f(t,y)
    depends on component j of y, and 0 otherwise.    
 
    ODE15S can solve problems M(t,y)*y' = f(t,y) with mass matrix M(t,y). Use
    ODESET to set the 'Mass' property to a function handle MASS if MASS(T,Y) 
    returns the value of the mass matrix. If the mass matrix is constant, 
    the matrix can be used as the value of the 'Mass' option. Problems with
    state-dependent mass matrices are more difficult. If the mass matrix does
    not depend on the state variable Y and the function MASS is to be called
    with one input argument T, set 'MStateDependence' to 'none'. If the mass
    matrix depends weakly on Y, set 'MStateDependence' to 'weak' (the
    default) and otherwise, to 'strong'. In either case the function MASS is
    to be called with the two arguments (T,Y). If there are many differential
    equations, it is important to exploit sparsity: Return a sparse
    M(t,y). Either supply the sparsity pattern of df/dy using the 'JPattern'
    property or a sparse df/dy using the Jacobian property. For strongly
    state-dependent M(t,y), set 'MvPattern' to a sparse matrix S with S(i,j)
    = 1 if for any k, the (i,k) component of M(t,y) depends on component j of
    y, and 0 otherwise.    
 
    If the mass matrix is non-singular, the solution of the problem is
    straightforward. See examples FEM1ODE, FEM2ODE, BATONODE, or
    BURGERSODE. If M(t0,y0) is singular, the problem is a differential-
    algebraic equation (DAE). ODE15S solves DAEs of index 1. DAEs have
    solutions only when y0 is consistent, i.e., there is a yp0 such that
    M(t0,y0)*yp0 = f(t0,y0). Use ODESET to set 'MassSingular' to 'yes', 'no',
    or 'maybe'. The default of 'maybe' causes ODE15S to test whether M(t0,y0)
    is singular. You can provide yp0 as the value of the 'InitialSlope'
    property. The default is the zero vector. If y0 and yp0 are not
    consistent, ODE15S treats them as guesses, tries to compute consistent
    values close to the guesses, and then goes on to solve the problem. See
    examples HB1DAE or AMP1DAE.  
 
    [TOUT,YOUT,TE,YE,IE] = ODE15S(ODEFUN,TSPAN,Y0,OPTIONS) with the 'Events'
    property in OPTIONS set to a function handle EVENTS, solves as above 
    while also finding where functions of (T,Y), called event functions, 
    are zero. For each function you specify whether the integration is 
    to terminate at a zero and whether the direction of the zero crossing 
    matters. These are the three column vectors returned by EVENTS: 
    [VALUE,ISTERMINAL,DIRECTION] = EVENTS(T,Y). For the I-th event function: 
    VALUE(I) is the value of the function, ISTERMINAL(I)=1 if the integration 
    is to terminate at a zero of this event function and 0 otherwise. 
    DIRECTION(I)=0 if all zeros are to be computed (the default), +1 if only 
    zeros where the event function is increasing, and -1 if only zeros where 
    the event function is decreasing. Output TE is a column vector of times 
    at which events occur. Rows of YE are the corresponding solutions, and 
    indices in vector IE specify which event occurred.    
    
    SOL = ODE15S(ODEFUN,[T0 TFINAL],Y0...) returns a structure that can be
    used with DEVAL to evaluate the solution or its first derivative at 
    any point between T0 and TFINAL. The steps chosen by ODE15S are returned 
    in a row vector SOL.x.  For each I, the column SOL.y(:,I) contains 
    the solution at SOL.x(I). If events were detected, SOL.xe is a row vector 
    of points at which events occurred. Columns of SOL.ye are the corresponding 
    solutions, and indices in vector SOL.ie specify which event occurred. 
 
    Example
          [t,y]=ode15s(@vdp1000,[0 3000],[2 0]);   
          plot(t,y(:,1));
      solves the system y' = vdp1000(t,y), using the default relative error
      tolerance 1e-3 and the default absolute tolerance of 1e-6 for each
      component, and plots the first component of the solution.
 
    See also ODE23S, ODE23T, ODE23TB, ODE45, ODE23, ODE113, ODE15I,
             ODESET, ODEPLOT, ODEPHAS2, ODEPHAS3, ODEPRINT, DEVAL,
             ODEEXAMPLES, VDPODE, BRUSSODE, HB1DAE, FUNCTION_HANDLE.

    Documentation for ode15s
       doc ode15s

help ode23t

 ODE23T  Solve moderately stiff ODEs and DAEs, trapezoidal rule.
    [TOUT,YOUT] = ODE23T(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates 
    the system of differential equations y' = f(t,y) from time T0 to TFINAL 
    with initial conditions Y0. ODEFUN is a function handle. For a scalar T 
    and a vector Y, ODEFUN(T,Y) must return a column vector corresponding 
    to f(t,y). Each row in the solution array YOUT corresponds to a time
    returned in the column vector TOUT.  To obtain solutions at specific
    times T0,T1,...,TFINAL (all increasing or all decreasing), use TSPAN =
    [T0 T1 ... TFINAL].     
    
    [TOUT,YOUT] = ODE23T(ODEFUN,TSPAN,Y0,OPTIONS) solves as above with default
    integration parameters replaced by values in OPTIONS, an argument created
    with the ODESET function. See ODESET for details. Commonly used options
    are scalar relative error tolerance 'RelTol' (1e-3 by default) and vector
    of absolute error tolerances 'AbsTol' (all components 1e-6 by default).  
    If certain components of the solution must be non-negative, use
    ODESET to set the 'NonNegative' property to the indices of these
    components.  The 'NonNegative' property is ignored for problems
    where there is a mass matrix. 
    
    The Jacobian matrix df/dy is critical to reliability and efficiency. Use
    ODESET to set 'Jacobian' to a function handle FJAC if FJAC(T,Y) returns 
    the Jacobian df/dy or to the matrix df/dy if the Jacobian is constant. 
    If the 'Jacobian' option is not set (the default), df/dy is approximated 
    by finite differences. Set 'Vectorized' 'on' if the ODE function is coded 
    so that ODEFUN(T,[Y1 Y2 ...]) returns [ODEFUN(T,Y1) ODEFUN(T,Y2) ...]. 
    If df/dy is a sparse matrix, set 'JPattern' to the sparsity pattern of
    df/dy, i.e., a sparse matrix S with S(i,j) = 1 if component i of f(t,y)
    depends on component j of y, and 0 otherwise.    
 
    ODE23T can solve problems M(t,y)*y' = f(t,y) with mass matrix M(t,y). Use
    ODESET to set the 'Mass' property to a function handle MASS if MASS(T,Y) 
    returns the value of the mass matrix. If the mass matrix is constant, 
    the matrix can be used as the value of the 'Mass' option. Problems with
    state-dependent mass matrices are more difficult. If the mass matrix does
    not depend on the state variable Y and the function MASS is to be called
    with one input argument T, set 'MStateDependence' to 'none'. If the mass
    matrix depends weakly on Y, set 'MStateDependence' to 'weak' (the
    default) and otherwise, to 'strong'. In either case the function MASS is
    to be called with the two arguments (T,Y). If there are many differential
    equations, it is important to exploit sparsity: Return a sparse
    M(t,y). Either supply the sparsity pattern of df/dy using the 'JPattern'
    property or a sparse df/dy using the Jacobian property. For strongly
    state-dependent M(t,y), set 'MvPattern' to a sparse matrix S with S(i,j)
    = 1 if for any k, the (i,k) component of M(t,y) depends on component j of
    y, and 0 otherwise.    
 
    If the mass matrix is non-singular, the solution of the problem is
    straightforward. See examples FEM1ODE, FEM2ODE, BATONODE, or
    BURGERSODE. If M(t0,y0) is singular, the problem is a differential-
    algebraic equation (DAE). ODE23T solves DAEs of index 1. DAEs have
    solutions only when y0 is consistent, i.e., there is a yp0 such that
    M(t0,y0)*yp0 = f(t0,y0). Use ODESET to set 'MassSingular' to 'yes', 'no',
    or 'maybe'. The default of 'maybe' causes ODE23T to test whether M(t0,y0)
    is singular. You can provide yp0 as the value of the 'InitialSlope'
    property. The default is the zero vector. If y0 and yp0 are not
    consistent, ODE23T treats them as guesses, tries to compute consistent
    values close to the guesses, and then goes on to solve the problem. See
    examples HB1DAE or AMP1DAE.    
 
    [TOUT,YOUT,TE,YE,IE] = ODE23T(ODEFUN,TSPAN,Y0,OPTIONS) with the 'Events'
    property in OPTIONS set to a function handle EVENTS, solves as above 
    while also finding where functions of (T,Y), called event functions, 
    are zero. For each function you specify whether the integration is 
    to terminate at a zero and whether the direction of the zero crossing 
    matters. These are the three column vectors returned by EVENTS: 
    [VALUE,ISTERMINAL,DIRECTION] = EVENTS(T,Y). For the I-th event function: 
    VALUE(I) is the value of the function, ISTERMINAL(I)=1 if the integration 
    is to terminate at a zero of this event function and 0 otherwise. 
    DIRECTION(I)=0 if all zeros are to be computed (the default), +1 if only 
    zeros where the event function is increasing, and -1 if only zeros where 
    the event function is decreasing. Output TE is a column vector of times 
    at which events occur. Rows of YE are the corresponding solutions, and 
    indices in vector IE specify which event occurred.    
 
    SOL = ODE23T(ODEFUN,[T0 TFINAL],Y0...) returns a structure that can be
    used with DEVAL to evaluate the solution or its first derivative at 
    any point between T0 and TFINAL. The steps chosen by ODE23T are returned 
    in a row vector SOL.x.  For each I, the column SOL.y(:,I) contains 
    the solution at SOL.x(I). If events were detected, SOL.xe is a row vector 
    of points at which events occurred. Columns of SOL.ye are the corresponding 
    solutions, and indices in vector SOL.ie specify which event occurred. 
 
    Example
          [t,y]=ode23t(@vdp1000,[0 3000],[2 0]);   
          plot(t,y(:,1));
      solves the system y' = vdp1000(t,y), using the default relative error
      tolerance 1e-3 and the default absolute tolerance of 1e-6 for each
      component, and plots the first component of the solution.
 
    See also ODE15S, ODE23S, ODE23TB, ODE45, ODE23, ODE113, ODE15I,
             ODESET, ODEPLOT, ODEPHAS2, ODEPHAS3, ODEPRINT, DEVAL,
             ODEEXAMPLES, VDPODE, BRUSSODE, HB1DAE, FUNCTION_HANDLE.

    Documentation for ode23t
       doc ode23t

If you want to fly, you need to use something designed to fly.