Hello, I have problem with odeToVectorField in ODE second order system of equations.

% clear all;
g= 9.81;
m= 0.468; 
Ix= 4.856e-3;
Iy= 4.856e-3;
Iz= 8.801e-3;
Ax= 0.25;
Ay= 0.25;
Az= 0.25;
l= 0.225;
k= 2.98e-6;
b= 1.14e-7;
Im= 3.357e-5;
f= 0; tu0= 0; tu1=0; tu2= 0; 
Tuo= [tu0 tu1 tu2];
syms x0(t) x1(t) x2(t) t0(t) t1(t) t2(t) t
X= [x0(t) x1(t) x2(t)];
T= [t0(t) t1(t) t2(t)];
dX= diff(X,t);
dT= diff(T,t);
% d2X= diff(dX,t);
% d2T= diff(dT,t);
s1= sin(t0) ;s2=sin(t1) ;s3=sin(t2);
c1= cos(t0) ;c2=cos(t1) ;c3=cos(t2);
t0d= dT(1); t1d=dT(2) ; t2d=dT(3) ;
% C Matrix
c11= 0;
c12= (Iy-Iz)*(t1d*c1*s1+ t2d*s1*s1*c2)+(Iz-Iy)*(t2d*c1*c1*c2)- Ix*t2d*c2;
c13= (Iz-Iy)*t2d*c1*s1*c2*c2;
c21= (Iz-Iy)*(t1d*c1*s1+t2d*s1*c2)+(Iy-Iz)*t2d*c1*c1*c2+Ix*t2d*c2;
c22= (Iz-Iy)*t0d*c1*s1 ;
c23= -Ix*t2d*s2*c2+Iy*t2d*s1*s1*s2*c2+Iz*t2d*c1*c1*s2*c2;
c31= (Iy-Iz)*t2d*c2*c2*s1*c1-Ix*t1d*c2;
c32= (Iz-Iy)*(t1d*c1*s1*s2+t0d*s1*s1*c2)+(Iy-Iz)*t0d*c1*c1*c2+...
    Ix*t2d*s2*c2-Iy*t2d*s1*s1*s2*c2-Iz*t2d*c1*c1*s2*c2;
c33= (Iy-Iz)*t0d*c1*s1*c2*c2-Iy*t1d*s1*s1*c2*s2-Iz*t1d*c1*c1*c2*s2+Ix*t1d*c2*s2;
C= [c11 c12 c13; c21 c22 c23; c31 c32 c33;];
% J matrix -invers
j11= Ix;
j12= 0;
j13= -Ix*s2;
j21= 0 ;
j22= Iy*c1*c1+Iz*s1*s1;
j23= (Iy-Iz)*c1*s1*c2;
j31= -Ix*s2;
j32= (Iy-Iz)*c1*s1*c2;
j33= Ix*s2*s2+Iy*s1*s1*c2*c2+Iz*c1*c1*c2*c2;
J = ([j11 j12 j13; j21 j22 j23; j31 j32 j33]);
DOM00= [c3*s2*c1+s3*s1, s3*s2*c1-c3*s1, c2*c1];
Ar= [Ax 0 0; 0 Ay 0; 0 0 Az];
%  Equation of Motion 
DOM01= J\(Tuo'- C*dT');
DOM02= [0 0 -g]' + (f/m)*DOM00'-(1.0/m)*Ar*dX';
Eq1 = diff(t0,2)== DOM01(1);
Eq2 = diff(t1,2)== DOM01(2);
Eq3 = diff(t2,2)== DOM01(3);
Eq4 = diff(x0,2)== DOM02(1);
Eq5 = diff(x1,2)== DOM02(2);
Eq6 = diff(x2,2)== DOM02(3);
vars = [t0(t); t1(t); t2(t); x0(t); x1(t); x2(t)];
[V,S] = odeToVectorField([Eq1 Eq2 Eq3 Eq4 Eq5 Eq6]);

Error using mupadengine/feval_internal
Unable to convert the initial value problem to an equivalent dynamical system. Either the differential equations cannot be solved for the highest derivatives or inappropriate initial conditions were specified.

Error in odeToVectorField>mupadOdeToVectorField (line 171)
T = feval_internal(symengine,'symobj::odeToVectorField',sys,x,stringInput);

Error in odeToVectorField (line 119)
sol = mupadOdeToVectorField(varargin);

M = matlabFunction(V,'vars', {'t'});

[Star Strider]

The ‘DOM01’ system appears to be nonlinear.

[Vahid]

thanks for reply

correct, it is not linear ...

[Roshana Sandadini]

Can I see the original Problem?

Answer 1

[Torsten]

dydt = fun(0.1,ones(12,1)).'

dydt = 1×12

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 3.7411 -0.6185 3.2759 -0.5342 -0.5342 -10.3442

function dydt = fun(t,y)
  g= 9.81;
  m= 0.468; 
  Ix= 4.856e-3;
  Iy= 4.856e-3;
  Iz= 8.801e-3;
  Ax= 0.25;
  Ay= 0.25;
  Az= 0.25;
  l= 0.225;
  k= 2.98e-6;
  b= 1.14e-7;
  Im= 3.357e-5;
  f= 0; tu0= 0; tu1=0; tu2= 0; 
  Tuo= [tu0 tu1 tu2];
  t0 = y(1);
  t1 = y(2);
  t2 = y(3);
  x0 = y(4);
  x1 = y(5);
  x2 = y(6);
  t0d = y(7);
  t1d = y(8);
  t2d = y(9);
  x0d = y(10);
  x1d = y(11);
  x2d = y(12);
  X = [x0 x1 x2];
  T = [t0 t1 t2];
  dX = [x0d x1d x2d];
  dT = [t0d t1d t2d];
  s1= sin(t0) ;s2=sin(t1) ;s3=sin(t2);
  c1= cos(t0) ;c2=cos(t1) ;c3=cos(t2);
  % C Matrix
  c11= 0;
  c12= (Iy-Iz)*(t1d*c1*s1+ t2d*s1*s1*c2)+(Iz-Iy)*(t2d*c1*c1*c2)- Ix*t2d*c2;
  c13= (Iz-Iy)*t2d*c1*s1*c2*c2;
  c21= (Iz-Iy)*(t1d*c1*s1+t2d*s1*c2)+(Iy-Iz)*t2d*c1*c1*c2+Ix*t2d*c2;
  c22= (Iz-Iy)*t0d*c1*s1 ;
  c23= -Ix*t2d*s2*c2+Iy*t2d*s1*s1*s2*c2+Iz*t2d*c1*c1*s2*c2;
  c31= (Iy-Iz)*t2d*c2*c2*s1*c1-Ix*t1d*c2;
  c32= (Iz-Iy)*(t1d*c1*s1*s2+t0d*s1*s1*c2)+(Iy-Iz)*t0d*c1*c1*c2+...
       Ix*t2d*s2*c2-Iy*t2d*s1*s1*s2*c2-Iz*t2d*c1*c1*s2*c2;
  c33= (Iy-Iz)*t0d*c1*s1*c2*c2-Iy*t1d*s1*s1*c2*s2-Iz*t1d*c1*c1*c2*s2+Ix*t1d*c2*s2;
  C= [c11 c12 c13; c21 c22 c23; c31 c32 c33;];
  % J matrix -invers
  j11= Ix;
  j12= 0;
  j13= -Ix*s2;
  j21= 0 ;
  j22= Iy*c1*c1+Iz*s1*s1;
  j23= (Iy-Iz)*c1*s1*c2;
  j31= -Ix*s2;
  j32= (Iy-Iz)*c1*s1*c2;
  j33= Ix*s2*s2+Iy*s1*s1*c2*c2+Iz*c1*c1*c2*c2;
  J = ([j11 j12 j13; j21 j22 j23; j31 j32 j33]);
  DOM00= [c3*s2*c1+s3*s1, s3*s2*c1-c3*s1, c2*c1];
  Ar= [Ax 0 0; 0 Ay 0; 0 0 Az];
  %  Equation of Motion 
  DOM01= J\(Tuo.'- C*dT.');
  DOM02= [0 0 -g].' + (f/m)*DOM00.'-(1.0/m)*Ar*dX.';
  dydt = [y(7:12);DOM01;DOM02];
end

Answer 2

[ProblemSolver]

The error you are getting indicates that the system of equations defined by DOM01 is non-linear. The odeToVectorField function in MATLAB requires the system of equations to be linear for it to work properly. But in your case, DOM01 contains trigonometric functions like sin() and cos() which makes it a nonlinear system.

Some options to resolve this:

1. Try linearizing the equations around an operating point using small angle approximations. This would eliminate the trig terms and result in a linear system that odeToVectorField can handle.
2. Use an alternate ODE solver in MATLAB that supports nonlinear systems, like ode45:

[T,X] = ode45(@ODEfunc,tspan,X0);

function dXdt = ODEfunc(t,X)
   % Implement RHS of equations
   dXdt = ... 
end

3. Convert the equations to state-space form with nonlinear dynamics:

   xdot = f(x,u)
   y = g(x,u)

And use Simulink to simulate the state-space model instead of solving analytically with odeToVectorField.

Answer 3

[Paul]

The problem is confusion between sym expressions and symfun objects.

g= 9.81;
m= 0.468; 
Ix= 4.856e-3;
Iy= 4.856e-3;
Iz= 8.801e-3;
Ax= 0.25;
Ay= 0.25;
Az= 0.25;
l= 0.225;
k= 2.98e-6;
b= 1.14e-7;
Im= 3.357e-5;
f= 0; tu0= 0; tu1=0; tu2= 0; 
Tuo= [tu0 tu1 tu2];
syms x0(t) x1(t) x2(t) t0(t) t1(t) t2(t) t
X= [x0(t) x1(t) x2(t)];
T= [t0(t) t1(t) t2(t)];
dX= diff(X,t);
dT= diff(T,t);
% d2X= diff(dX,t);
% d2T= diff(dT,t);
s1= sin(t0) ;s2=sin(t1) ;s3=sin(t2);
c1= cos(t0) ;c2=cos(t1) ;c3=cos(t2);
t0d= dT(1); t1d=dT(2) ; t2d=dT(3) ;
% C Matrix
c11= 0;
c12= (Iy-Iz)*(t1d*c1*s1+ t2d*s1*s1*c2)+(Iz-Iy)*(t2d*c1*c1*c2)- Ix*t2d*c2;
c13= (Iz-Iy)*t2d*c1*s1*c2*c2;
c21= (Iz-Iy)*(t1d*c1*s1+t2d*s1*c2)+(Iy-Iz)*t2d*c1*c1*c2+Ix*t2d*c2;
c22= (Iz-Iy)*t0d*c1*s1 ;
c23= -Ix*t2d*s2*c2+Iy*t2d*s1*s1*s2*c2+Iz*t2d*c1*c1*s2*c2;
c31= (Iy-Iz)*t2d*c2*c2*s1*c1-Ix*t1d*c2;
c32= (Iz-Iy)*(t1d*c1*s1*s2+t0d*s1*s1*c2)+(Iy-Iz)*t0d*c1*c1*c2+...
    Ix*t2d*s2*c2-Iy*t2d*s1*s1*s2*c2-Iz*t2d*c1*c1*s2*c2;
c33= (Iy-Iz)*t0d*c1*s1*c2*c2-Iy*t1d*s1*s1*c2*s2-Iz*t1d*c1*c1*c2*s2+Ix*t1d*c2*s2;
C= [c11 c12 c13; c21 c22 c23; c31 c32 c33;];
% J matrix -invers
j11= Ix;
j12= 0;
j13= -Ix*s2;
j21= 0 ;
j22= Iy*c1*c1+Iz*s1*s1;
j23= (Iy-Iz)*c1*s1*c2;
j31= -Ix*s2;
j32= (Iy-Iz)*c1*s1*c2;
j33= Ix*s2*s2+Iy*s1*s1*c2*c2+Iz*c1*c1*c2*c2;
J = ([j11 j12 j13; j21 j22 j23; j31 j32 j33]);
DOM00= [c3*s2*c1+s3*s1, s3*s2*c1-c3*s1, c2*c1];
Ar= [Ax 0 0; 0 Ay 0; 0 0 Az];
%  Equation of Motion 

Using ctranspose here will result in a lot of unnecessary conj() calls later on. Use transpose.

%DOM01= J\(Tuo'- C*dT');
%DOM02= [0 0 -g]' + (f/m)*DOM00'-(1.0/m)*Ar*dX';
DOM01= J\(Tuo.'- C*dT.');
DOM02= [0 0 -g].' + (f/m)*DOM00.'-(1.0/m)*Ar*dX.';

At this point, DOM01 is scalar symfun that takes in a scalar value and returns a 3 x 1 sym expression

whos DOM01

  Name       Size            Bytes  Class     Attributes

  DOM01      1x1                 8  symfun              

temp = DOM01(0);
whos temp

  Name      Size            Bytes  Class    Attributes

  temp      3x1                 8  sym                

So the following assignment is evaluating DOM01 at t = 1 and returning the 3 x 1 result and setting that equal to the 3x1 expansion of the lhs. The rhs is too painful to look at, but if you do you'll see various terms being evaluated at t = 1.

Eq1 = diff(t0,2)== DOM01(1);
lhs(Eq1)

ans(t) = 

I suspect that you really want to pick off the elements of DOM01, like so

Eq1 = diff(t0,2) == [1 0 0]*DOM01(t)

Eq1(t) = 

You could do that for all of the equations. Or just do it all at once:

%{
Eq1 = diff(t0,2)== DOM01(1);
Eq2 = diff(t1,2)== DOM01(2);
Eq3 = diff(t2,2)== DOM01(3);
Eq4 = diff(x0,2)== DOM02(1);
Eq5 = diff(x1,2)== DOM02(2);
Eq6 = diff(x2,2)== DOM02(3);
vars = [t0(t); t1(t); t2(t); x0(t); x1(t); x2(t)];
[V,S] = odeToVectorField([Eq1 Eq2 Eq3 Eq4 Eq5 Eq6]);
%}
Eq = [diff(t0,2);diff(t1,2);diff(t2,2);diff(x0,2);diff(x1,2);diff(x2,2)] == [DOM01(t);DOM02(t)];
[V,S] = odeToVectorField(Eq);  % no errors

odefun = matlabFunction(V,'vars',{'t','Y'})

odefun = function_handle with value:
    @(t,Y)[Y(2);((cos(Y(3)).*sin(Y(3)).*Y(2).*Y(4).*-2.546391437701197e+31-cos(Y(1)).*cos(Y(3)).^2.*sin(Y(1)).*Y(6).^2.*5.680808882942518e+31+cos(Y(1)).*cos(Y(3)).^4.*sin(Y(1)).*Y(6).^2.*1.029588117355377e+32-cos(Y(1)).*sin(Y(1)).*sin(Y(3)).^2.*Y(6).^2.*3.134417445241321e+31+cos(Y(1)).*sin(Y(1)).*sin(Y(3)).^4.*Y(6).^2.*3.134417445241321e+31+cos(Y(1)).*cos(Y(3)).^2.*Y(4).*Y(6).*5.680808882942518e+31-cos(Y(1)).*cos(Y(3)).^4.*Y(4).*Y(6).*4.615072290611251e+31+cos(Y(3)).*sin(Y(3)).^3.*Y(2).*Y(4).*7.161463728312448e+31+cos(Y(3)).^3.*sin(Y(3)).*Y(2).*Y(4).*7.161463728312448e+31+cos(Y(1)).*sin(Y(3)).^2.*Y(4).*Y(6).*3.134417445241321e+31+cos(Y(1)).*sin(Y(3)).^3.*Y(4).*Y(6).*2.546391437701197e+31+cos(Y(3)).*sin(Y(1)).*sin(Y(3)).*Y(2).*Y(6).*2.546391437701197e+31+cos(Y(1)).*cos(Y(3)).^2.*sin(Y(3)).*Y(4).*Y(6).*4.615072290611251e+31-cos(Y(3)).*sin(Y(1)).*sin(Y(3)).^3.*Y(2).*Y(6).*7.161463728312448e+31-cos(Y(3)).^3.*sin(Y(1)).*sin(Y(3)).*Y(2).*Y(6).*7.161463728312448e+31+cos(Y(1)).*cos(Y(3)).^2.*sin(Y(1)).*sin(Y(3)).^2.*Y(6).^2.*1.343029861879509e+32-cos(Y(1)).*cos(Y(3)).^2.*sin(Y(3)).^2.*Y(4).*Y(6).*6.683753143521304e+31).*(-1.760312695965974e-32))./(cos(Y(3)).^2.*sin(Y(3)).^2.*2.0+cos(Y(3)).^4+sin(Y(3)).^4);Y(4);(cos(Y(1)).*cos(Y(3)).*sin(Y(3)).^5.*Y(2).^2.*4.615072290611251e+31+cos(Y(1)).*cos(Y(3)).^5.*sin(Y(3)).*Y(2).^2.*4.615072290611251e+31+cos(Y(3)).^2.*sin(Y(1)).*Y(2).*Y(4).*3.134417445241321e+31+cos(Y(3)).^4.*sin(Y(1)).*Y(2).*Y(4).*2.546391437701197e+31+sin(Y(1)).*sin(Y(3)).^2.*Y(2).*Y(4).*5.680808882942518e+31-sin(Y(1)).*sin(Y(3)).^4.*Y(2).*Y(4).*4.615072290611251e+31+cos(Y(1)).*cos(Y(3)).^3.*sin(Y(3)).^3.*Y(2).^2.*9.230144581222501e+31-cos(Y(1)).^3.*cos(Y(3)).*sin(Y(3)).^5.*Y(6).^2.*4.615072290611251e+31-cos(Y(1)).^3.*cos(Y(3)).^5.*sin(Y(3)).*Y(6).^2.*4.615072290611251e+31+cos(Y(1)).^2.*cos(Y(3)).^4.*Y(2).*Y(6).*5.680808882942518e+31-cos(Y(1)).^2.*cos(Y(3)).^6.*Y(2).*Y(6).*4.615072290611251e+31-cos(Y(3)).^2.*sin(Y(1)).^2.*Y(2).*Y(6).*3.134417445241321e+31+cos(Y(1)).^2.*sin(Y(3)).^4.*Y(2).*Y(6).*5.680808882942518e+31+cos(Y(3)).^4.*sin(Y(1)).^2.*Y(2).*Y(6).*8.81522632818384e+31+cos(Y(1)).^2.*sin(Y(3)).^6.*Y(2).*Y(6).*4.615072290611251e+31-sin(Y(1)).^2.*sin(Y(3)).^2.*Y(2).*Y(6).*5.680808882942518e+31+sin(Y(1)).^2.*sin(Y(3)).^4.*Y(2).*Y(6).*1.597669005649629e+32-cos(Y(1)).^3.*cos(Y(3)).^3.*sin(Y(3)).^3.*Y(6).^2.*9.230144581222501e+31-cos(Y(3)).^2.*sin(Y(1)).*sin(Y(3)).^2.*Y(2).*Y(4).*2.068680852910054e+31+cos(Y(1)).^2.*cos(Y(3)).^2.*sin(Y(3)).^2.*Y(2).*Y(6).*1.136161776588504e+32+cos(Y(1)).^2.*cos(Y(3)).^2.*sin(Y(3)).^4.*Y(2).*Y(6).*4.615072290611251e+31-cos(Y(1)).^2.*cos(Y(3)).^4.*sin(Y(3)).^2.*Y(2).*Y(6).*4.615072290611251e+31+cos(Y(3)).^2.*sin(Y(1)).^2.*sin(Y(3)).^2.*Y(2).*Y(6).*2.479191638468013e+32+cos(Y(1)).*cos(Y(3)).*sin(Y(1)).^2.*sin(Y(3)).*Y(6).^2.*2.546391437701197e+31-cos(Y(1)).*cos(Y(3)).*sin(Y(1)).^2.*sin(Y(3)).^3.*Y(6).^2.*7.161463728312448e+31-cos(Y(1)).*cos(Y(3)).^3.*sin(Y(1)).^2.*sin(Y(3)).*Y(6).^2.*7.161463728312448e+31-cos(Y(1)).*cos(Y(3)).*sin(Y(1)).*sin(Y(3)).*Y(4).*Y(6).*2.546391437701197e+31-cos(Y(1)).*cos(Y(3)).*sin(Y(1)).*sin(Y(3)).^2.*Y(4).*Y(6).*2.068680852910054e+31+cos(Y(1)).*cos(Y(3)).*sin(Y(1)).*sin(Y(3)).^3.*Y(4).*Y(6).*9.230144581222501e+31+cos(Y(1)).*cos(Y(3)).^3.*sin(Y(1)).*sin(Y(3)).*Y(4).*Y(6).*7.161463728312448e+31)./(cos(Y(1)).*(cos(Y(3)).^2.*sin(Y(3)).^2.*2.0+cos(Y(3)).^4+sin(Y(3)).^4).*5.680808882942518e+31);Y(6);(cos(Y(3)).^2.*Y(2).*Y(4).*3.134417445241321e+31+cos(Y(3)).^4.*Y(2).*Y(4).*2.546391437701197e+31+sin(Y(3)).^2.*Y(2).*Y(4).*5.680808882942518e+31-sin(Y(3)).^4.*Y(2).*Y(4).*4.615072290611251e+31-cos(Y(3)).^2.*sin(Y(1)).*Y(2).*Y(6).*3.134417445241321e+31+cos(Y(3)).^4.*sin(Y(1)).*Y(2).*Y(6).*8.81522632818384e+31-sin(Y(1)).*sin(Y(3)).^2.*Y(2).*Y(6).*5.680808882942518e+31+sin(Y(1)).*sin(Y(3)).^4.*Y(2).*Y(6).*1.597669005649629e+32-cos(Y(3)).^2.*sin(Y(3)).^2.*Y(2).*Y(4).*2.068680852910054e+31+cos(Y(3)).^2.*sin(Y(1)).*sin(Y(3)).^2.*Y(2).*Y(6).*2.479191638468013e+32+cos(Y(1)).*cos(Y(3)).*sin(Y(1)).*sin(Y(3)).*Y(6).^2.*2.546391437701197e+31-cos(Y(1)).*cos(Y(3)).*sin(Y(3)).*Y(4).*Y(6).*2.546391437701197e+31-cos(Y(1)).*cos(Y(3)).*sin(Y(1)).*sin(Y(3)).^3.*Y(6).^2.*7.161463728312448e+31-cos(Y(1)).*cos(Y(3)).^3.*sin(Y(1)).*sin(Y(3)).*Y(6).^2.*7.161463728312448e+31-cos(Y(1)).*cos(Y(3)).*sin(Y(3)).^2.*Y(4).*Y(6).*2.068680852910054e+31+cos(Y(1)).*cos(Y(3)).*sin(Y(3)).^3.*Y(4).*Y(6).*9.230144581222501e+31+cos(Y(1)).*cos(Y(3)).^3.*sin(Y(3)).*Y(4).*Y(6).*7.161463728312448e+31)./(cos(Y(1)).*(cos(Y(3)).^2.*sin(Y(3)).^2.*2.0+cos(Y(3)).^4+sin(Y(3)).^4).*5.680808882942518e+31);Y(8);Y(8).*(-1.25e+2./2.34e+2);Y(10);Y(10).*(-1.25e+2./2.34e+2);Y(12);Y(12).*(-1.25e+2./2.34e+2)-9.81e+2./1.0e+2]

Also, you may want consider defining all those constants as sym constants up front, like

g = sym(9.81);

In general, it's best (IMO) to stay in sym world as long as possible to allow the engine to find cancellations and simplifications, should there be any.

[Walter Roberson]

If you are going to be using floating point constants in symbolic contexts, then I recommend converting them to symbolic explicitly, so that you have control over the symbolic representation.

There is an example at https://www.mathworks.com/matlabcentral/answers/2009892-errors-in-using-null-command-due-to-truncation-error#comment_2853442 of a case where converting the numbers before using them in a symbolic expression was necessary in order for the mathematics of the situation to work out.

[Torsten]

@Paul

Why are

t0d= dT(1); t1d=dT(2) ; t2d=dT(3)

not 1x3 sym each representing dT evaluated at t = 1, t = 2 and t = 3 ?

[Paul]

g= 9.81;
m= 0.468; 
Ix= 4.856e-3;
Iy= 4.856e-3;
Iz= 8.801e-3;
Ax= 0.25;
Ay= 0.25;
Az= 0.25;
l= 0.225;
k= 2.98e-6;
b= 1.14e-7;
Im= 3.357e-5;
f= 0; tu0= 0; tu1=0; tu2= 0; 
Tuo= [tu0 tu1 tu2];
syms x0(t) x1(t) x2(t) t0(t) t1(t) t2(t) t
X= [x0(t) x1(t) x2(t)];
T= [t0(t) t1(t) t2(t)];
dX= diff(X,t);
dT= diff(T,t);
% d2X= diff(dX,t);
% d2T= diff(dT,t);
s1= sin(t0) ;s2=sin(t1) ;s3=sin(t2);
c1= cos(t0) ;c2=cos(t1) ;c3=cos(t2);
t0d= dT(1); t1d=dT(2) ; t2d=dT(3) ;

dT is a 1 x 3 sym, not a symfun

whos dT

  Name      Size            Bytes  Class    Attributes

  dT        1x3                 8  sym                

So dT(ii) is ordinary, vector indexing and returns what we want

dT(2)  % for example

ans = 

It would be a different story had it been coded like this

dT(t) = diff(T,t)  % dT is now a symfun with one argument that returns a 1 x 3 sym

dT(t) = 

whos dT

  Name      Size            Bytes  Class     Attributes

  dT        1x1                 8  symfun              

dT(1)  % oops!

ans = 

whos ans

  Name      Size            Bytes  Class    Attributes

  ans       1x3                 8  sym                

Surprisingly, diff does not have the standard 'Output Arguments' section to explain how this works. Now that I'm poking around, maybe it's not so surprising. int also does not have that section.

[Torsten]

@Paul

That's really subtle. Thank you for the explanation.

[Paul]

@Walter Roberson

To your point, is it possible that using symbolic constants might not always work. That is if we do this

syms a b
temp = a*b;
c = temp/a

c = 

b

the 'a's will always cancel.

If a and b are assigned symbolic constants

a = sym(1.23456);
b = sym(5.43210);
temp = a*b

temp = 

c = temp/a

c = 

then in this case c is returned with the value of b. Will this result always obtain because the symbolic constants are always rational numbers and so cancelling the common factors in computations is straightforward? Or can there be a case, say where c is result of a sequence of expressions involving who knows what kind of operations, where something gets lost in the shuffle, so to speak, and the result for c would be not quite what would be expected?

[Walter Roberson]

sym() of a floating point number can return multiples of π and square roots when the default conversion is used. The default conversion also approximates... so for example a floating point number sufficiently close to 1/3 will be converted to 1/3, not just the single double precision number that most closely represents 1/3. π ± about 999 eps is converted to π

But once you have converted to rational or sqrt or pi, the symbolic engine sticks with the representation and calculates with it until you have a double() or vpa() or vpaintegral(), or you have a relational operation in a context that forces a decision to be made. "if A<B" can be difficult to answer when A and B are symbolic formula. The engine cannot rely on being able to come up with mathematical proofs so sometimes it has to resort to numeric evaluation. But it holds off until then.

For example sum(sin(sym(linspace(0,pi))) will give a mostly unevaluated sum of hundreds of sin terms: it does not approximate the sin of rationals.

Now if you construct a symbolic floating-point number then the engine will create approximations as it goes. for example sym('1.23') or str2sym('1.23') or sym() of a floating point number with particular options to the call. Once you are working with symbolic floating-point including if you vpa() then you can run into the same kind of issues as you can with double precision