Error while using DAE example

I'm trying out the following example, http://in.mathworks.com/help/symbolic/set-up-your-dae-problem.html, of simulating a pendulum, whose characteristics is described by a differential algebraic expression. But I've hit a dead end at step 5, finding consistent initial consitions. For some reason the decic function returns the following error. I'm using matlab 2014b.

Error using
symengine>makeFhandle/@(t,in2,in3,param1,param2,param3)[-(in2(3,:).*in2(1,:)-in2(7,:).*param1.*param2)./param2;(-in2(3,:).*in2(2,:)+in2(6,:).*param1.*param2+param1.*param2.*param3)./param2;-param2.^2+in2(1,:).^2+in2(2,:).^2;in2(4,:).*in2(1,:).*2.0+in2(5,:).*in2(2,:).*2.0;in2(7,:).*in2(1,:).*2.0+in3(5,:).*in2(2,:).*2.0+in2(4,:).^2.*2.0+in2(5,:).^2.*2.0;in2(6,:)-in3(5,:);in2(5,:)-in3(2,:)]
Not enough input arguments.

Error in decic (line 66)
res = feval(odefun,t0,y0,yp0,varargin{:});

Error in Untitled (line 38)
[y0, yp0] = decic(f, 0, y0est, [], yp0est, [], opt)

Here's the entire code I tried.

clear all
clc

%Equations and variables
syms x(t) y(t) T(t) m r g;
eqs= [m*diff(x(t), 2) == T(t)/r*x(t), ...
      m*diff(y(t), 2) == T(t)/r*y(t) - m*g, ...
      x(t)^2 + y(t)^2 == r^2];

vars = [x(t); y(t); T(t)];

%reducing differential order

[eqs, vars, R] = reduceDifferentialOrder(eqs, vars);

%reduce DAE index
[DAEs,DAEvars] = reduceDAEIndex(eqs,vars);
[DAEs,DAEvars] = reduceRedundancies(DAEs,DAEvars)

%Finding function handles such as F(t,y(t),y'(t))=0 is reqd for ode15i, So
%convert a DAE system of function handle F=F(t,y,yp), where t is scalar , y
% yp vectors.
f = daeFunction(DAEs, DAEvars, m, r, g);
m = 1.0;
r = 1.0;
g = 9.81;
%The function handle f still contains symbolic parameters. Create a purely
%numeric function handle F that you can pass to ode15i.
F = @(t, Y, YP) f(t, Y, YP, m, r, g);

DAEs = subs(DAEs);
%F = daeFunction(DAEs, DAEvars);

%finding consistent initial conditions
y0est = [0.5*r; -0.8*r; 0; 0; 0; 0; 0];
yp0est = zeros(7,1);
opt = odeset('RelTol', 10.0^(-7), 'AbsTol' , 10.0^(-7));
[y0, yp0] = decic(f, 0, y0est, [], yp0est, [], opt)

Answer 1

[Walter Roberson]

The documentation for decic has no explicit information about how many parameters the odefun takes. But it refers to the documentation for ode15i and that has an example that shows the undefined function @weissinger being invoked by both decic and ode15i, so the odefun for decic must that the same parameters as the odefun for ode15i . Which the ode15i documentation does not specify. But the documentation talking about the ode solvers says, in the description of the basic input arguments, that the form for odefun is dydt = odefun(t,y) . Therefore odefun must take two input parameters, and therefore decic must be called on a function with two input parameters. The counter-argument is that the place of the crash shows

res = feval(odefun,t0,y0,yp0,varargin{:});

which implies odefun is being called with 3 arguments.

You are calling decic on "f", which you created via

f = daeFunction(DAEs, DAEvars, m, r, g);

but notice the comments after that indicating that

%The function handle f still contains symbolic parameters. Create a purely
%numeric function handle F that you can pass to ode15i.

Hence when you are passing f to decic you are passing a symbolic expression. feval() is trying to convert that symbolic expression into a function handle, but there are 6 symbols in the symbolic expression so the function handle is being created with 6 parameters. But decic is only passing 3 parameters to it.

It appears to me that you should be passing F to decic instead of f.

[HJ Jay]

Thanks, Walter. You're right.