Error using feval Function to evaluate must be represented as a string scalar, character vector, or function_handle object.

% Example5_5.m
% Solution to Example 5.5. This program calculates and plots
% the concentration of cell mass, concentration of penicillin,
% optimal temperature profile, and adjoint variable of a batch
% penicillin fermentor. It uses the function COLLOCATION to
% solve the set of system and adjoint equations.
clear
clc
clf
% Input data
w = input(‘ Enter w”s as a vector : ‘);
y0 = input(‘ Vector of known initial conditions = ‘);
yf = input(‘ Vector of final conditions = ‘);
guess = input(‘ Vector of guessed initial conditions = ‘);
fname = input(‘n M-file containing the set of differential equations : ‘);
fth = input(‘ M-file containing the necessary condition function : ‘);
n = input(‘ Number of internal collocation points = ‘);
rho = input(‘ Relaxation factor = ‘);
% Solution of the set of differential equations
[t,y] = collocation(fname,0,1,y0,yf,guess,n,rho,[],w,fth);
% Temperature changes
for k = 1:n+2
options=optimset;
theta(k) = fzero(fth,30,options,y(:,k),w);
end
% Plotting the results
subplot(2,2,1), plot(t,y(1,:))
xlabel(‘Time’)
ylabel(‘Cell’)
title(‘(a)’)
subplot(2,2,2), plot(t,y(2,:))
xlabel(‘Time’)
ylabel(‘Penicillin’)
title(‘(b)’)
subplot(2,2,3), plot(t,y(3,:))
xlabel(‘Time’)
ylabel(‘First Adjoint’)
title(‘(c)’)
subplot(2,2,4), plot(t,theta)
xlabel(‘Time’)
ylabel(‘Temperature (deg C)’)
title(‘(d)’)

function f = Ex5_5_func(t,y,w,fth)
% Function Ex5_5_func.M
% This function introduces the set of ordinary differential
% equations used in Example 5.5.
% Temperature
options=optimset;
theta = fzero(fth,30,options,y,w);
% Calculating the b’s
b1 = w(1) * (1-w(2)*(theta-w(3))^2) / (1-w(2)*(25-w(3))^2);
if b1<0, b1=0; end
b2 = w(4) * (1-w(2)*(theta-w(3))^2) / (1-w(2)*(25-w(3))^2);
if b2<0, b2=1e-6; end
b3 = w(5) * (1-w(2)*(theta-w(6))^2) / (1-w(2)*(25-w(6))^2);
if b3<0, b3=0; end
% Evaluating the function values
f(1) = b1*y(1) – b1/b2*y(1)^2;
f(2) = b3*y(1);
f(3) = -b1*y(3) + 2*b1/b2*y(1)*y(3) – b3;
f = f’; % Make it a column vector

function ftheta = Ex5_5_theta(theta,y,w)
% Function Ex5_5_theta.M
% This function calculates the value of the necessary condition
% as a function of the temperature (theta). It is used in solving
% Example 5.5.
% Calculating the b’s
b1 = w(1) * (1-w(2)*(theta-w(3))^2) / (1-w(2)*(25-w(3))^2);
db1 = w(1)*(-w(2))*2*(theta-w(3)) / (1-w(2)*(25-w(3))^2);
b2 = w(4) * (1-w(2)*(theta-w(3))^2) / (1-w(2)*(25-w(3))^2);
db2 = w(4)*(-w(2))*2*(theta-w(3)) / (1-w(2)*(25-w(3))^2);
b3 = w(5) * (1-w(2)*(theta-w(6))^2) / (1-w(2)*(25-w(6))^2);
db3 = w(5)*(-w(2))*2*(theta-w(6)) / (1-w(2)*(25-w(6))^2);
% The function
ftheta = y(3)*(y(1)*db1-y(1)^2*(db1*b2-db2*b1)/b2^2)+y(1)*db3;

function [x,y] = collocation(ODEfile,x0,xf,y0,yf,guess,n,rho,tol,varargin)
%COLLOCATION Solves a boundary value set of ordinary differential
% equations by the orthogonal collocation method.
%
% [X,Y]=COLLOCATION(‘F’,X0,XF,Y0,YF,GAMMA,N) integrates the set of
% ordinary differential equations from X0 to XF by the Nth-degree
% orthogonal collocation method. The equations are contained in
% the M-file F.M. Y0, YF, and GAMMA are the vectors of initial
% conditions, final conditions, and starting guesses respectively.
% The function returns the independent variable in the vector X
% and the set of dependent variables in the matrix Y.
%
% [X,Y]=COLLOCATION(‘F’,X0,XF,Y0,YF,GAMMA,N,RHO,TOL,P1,P2,…)
% uses relaxation factor RHO and tolerance TOL for convergence
% test. Additional parameters P1, P2, … are passed directly to
% the function F. Pass an empty matrix for RHO or TOL to use the
% default value.
%
% See also SHOOTING
% (c) N. Mostoufi & A. Constantinides
% January 1, 1999
% Initialization
if nargin < 7 | isempty(n)
n = 1;
end
if nargin < 8 | isempty(rho)
rho = 1;
end
if nargin < 9 | isempty(tol)
tol = 1e-6;
end
y0 = (y0(:).’)’; % Make sure it’s a column vector
yf = (yf(:).’)’; % Make sure it’s a column vector
guess = (guess(:).’)’; % Make sure it’s a column vector
% Checking the number of guesses
if length(yf) ~= length(guess)
error(‘ The number of guessed conditions is not equal to the number of final conditions.’)
end
r = length(y0); % Number of initial conditions
m = r + length(yf); % Number of boundary conditions
% Checking the number of equations
ftest = feval(ODEfile,x0,[y0 ; guess],varargin{:});
if length(ftest) ~= m
error(‘ The number of equations is not equal to the number of boundary conditions.’)
end
fprintf(‘n Integrating. Please wait.nn’)
% Coefficients of the Legendre polynomial
for k = 0 : n/2
cl(2*k+1) = (-1)^k * gamma(2*n-2*k+1) / …
(2^n * gamma(k+1) * gamma(n-k+1) * gamma(n-2*k+1));
if k < n/2
cl(2*k+2) = 0;
end
end
zl = roots(cl); % Roots of the Legendre polynomial
z = [-1; sort(zl); 1]; % Collocation points (z)
x = (xf-x0)*z/2+(xf+x0)/2; % Collocation points (x)
% Bulding the vector of starting values of the dependent variables
[p,q] = RK(ODEfile,x0,xf,(xf-x0)/20,[y0 ; guess],2,varargin{:});
for k = 1:m
y(k,:) = spline(p,q(k,:),x’);
end
y(r+1:m,end) = yf(1:m-r);
% Building the matrix A
Q(:,1) = ones(n+2,1);
C(:,1) = zeros(n+2,1);
for i = 1:n+1
Q(:,i+1) = x.^i;
C(:,i+1) = i*x.^(i-1);
end
A = C*inv(Q);
for k = 1:m
k1 = (k-1)*(n+2)+1;
k2 = k1 + n+1;
Am(k1:k2,k1:k2) = A; % Building the matrix Am
Y(k1:k2) = y(k,:); % Building the vector Y
end
Y = Y’; % Make it a column vector
Y1 = Y * 1.1;
iter = 0;
maxiter = 100;
F = zeros(m*(n+2),1);
Fa = zeros(m*(n+2),1);
dY = zeros(m*(n+2),1);
position = []; % Collocation points excluding boundary conditions
for k = 1:m
if k <= r
position = [position, (k-1)*(n+2)+[2:n+2] ];
else
position = [position, (k-1)*(n+2)+[1:n+1] ];
end
end
% Newton’s method
while max(abs(Y1 – Y)) > tol & iter < maxiter
iter = iter + 1;
fprintf(‘ Iteration %3dn’,iter)
Y1 = Y;
% Building the vector F
for k = 1:n+2
F(k : n+2 : (m-1)*(n+2)+k) = …
feval(ODEfile,x(k),Y(k : n+2 : (m-1)*(n+2)+k),varargin{:});
end
fnk = Am * Y – F;

% Set dY for derivation
for k = 1:m*(n+1)
if Y(position(k)) ~= 0
dY(position(k)) = Y(position(k)) / 100;
else
dY(position(k)) = 0.01;
end
end

% Calculation of the Jacobian matrix
for k = 1:m
for kk = 1:n+1
a = Y;
nc = (k-1)*(n+1)+kk;
a(position(nc)) = Y(position(nc)) + dY(position(nc));
for kkk = 1:n+2
Fa(kkk : n+2 : (m-1)*(n+2)+kkk) = …
feval(ODEfile,x(kkk),a(kkk:n+2:(m-1)*(n+2)+kkk),varargin{:});
end
fnka = Am * a – Fa;
jacob(:,nc) = (fnka(position) – fnk(position))…
/ dY(position(nc));
end
end
% Next approximation of the roots
if det(jacob) == 0
Y(position) = Y(position) + max([abs(dY(position)); 1.1*tol]);
else
Y(position) = Y(position) – rho * inv(jacob) * fnk(position);
end
end
% Rearranging the y’s
for k = 1:m
k1 = (k-1)*(n+2)+1;
k2 = k1 + n+1;
y(k,:) = Y(k1:k2)’;
end
x = x’;
if iter >= maxiter
disp(‘Warning : Maximum iterations reached.’)
end% Example5_5.m
% Solution to Example 5.5. This program calculates and plots
% the concentration of cell mass, concentration of penicillin,
% optimal temperature profile, and adjoint variable of a batch
% penicillin fermentor. It uses the function COLLOCATION to
% solve the set of system and adjoint equations.
clear
clc
clf
% Input data
w = input(‘ Enter w”s as a vector : ‘);
y0 = input(‘ Vector of known initial conditions = ‘);
yf = input(‘ Vector of final conditions = ‘);
guess = input(‘ Vector of guessed initial conditions = ‘);
fname = input(‘n M-file containing the set of differential equations : ‘);
fth = input(‘ M-file containing the necessary condition function : ‘);
n = input(‘ Number of internal collocation points = ‘);
rho = input(‘ Relaxation factor = ‘);
% Solution of the set of differential equations
[t,y] = collocation(fname,0,1,y0,yf,guess,n,rho,[],w,fth);
% Temperature changes
for k = 1:n+2
options=optimset;
theta(k) = fzero(fth,30,options,y(:,k),w);
end
% Plotting the results
subplot(2,2,1), plot(t,y(1,:))
xlabel(‘Time’)
ylabel(‘Cell’)
title(‘(a)’)
subplot(2,2,2), plot(t,y(2,:))
xlabel(‘Time’)
ylabel(‘Penicillin’)
title(‘(b)’)
subplot(2,2,3), plot(t,y(3,:))
xlabel(‘Time’)
ylabel(‘First Adjoint’)
title(‘(c)’)
subplot(2,2,4), plot(t,theta)
xlabel(‘Time’)
ylabel(‘Temperature (deg C)’)
title(‘(d)’)

% Set dY for derivation
for k = 1:m*(n+1)
if Y(position(k)) ~= 0
dY(position(k)) = Y(position(k)) / 100;
else
dY(position(k)) = 0.01;
end
end

% Calculation of the Jacobian matrix
for k = 1:m
for kk = 1:n+1
a = Y;
nc = (k-1)*(n+1)+kk;
a(position(nc)) = Y(position(nc)) + dY(position(nc));
for kkk = 1:n+2
Fa(kkk : n+2 : (m-1)*(n+2)+kkk) = …
feval(ODEfile,x(kkk),a(kkk:n+2:(m-1)*(n+2)+kkk),varargin{:});
end
fnka = Am * a – Fa;
jacob(:,nc) = (fnka(position) – fnk(position))…
/ dY(position(nc));
end
end
% Next approximation of the roots
if det(jacob) == 0
Y(position) = Y(position) + max([abs(dY(position)); 1.1*tol]);
else
Y(position) = Y(position) – rho * inv(jacob) * fnk(position);
end
end
% Rearranging the y’s
for k = 1:m
k1 = (k-1)*(n+2)+1;
k2 = k1 + n+1;
y(k,:) = Y(k1:k2)’;
end
x = x’;
if iter >= maxiter
disp(‘Warning : Maximum iterations reached.’)
end % Example5_5.m
% Solution to Example 5.5. This program calculates and plots
% the concentration of cell mass, concentration of penicillin,
% optimal temperature profile, and adjoint variable of a batch
% penicillin fermentor. It uses the function COLLOCATION to
% solve the set of system and adjoint equations.
clear
clc
clf
% Input data
w = input(‘ Enter w”s as a vector : ‘);
y0 = input(‘ Vector of known initial conditions = ‘);
yf = input(‘ Vector of final conditions = ‘);
guess = input(‘ Vector of guessed initial conditions = ‘);
fname = input(‘n M-file containing the set of differential equations : ‘);
fth = input(‘ M-file containing the necessary condition function : ‘);
n = input(‘ Number of internal collocation points = ‘);
rho = input(‘ Relaxation factor = ‘);
% Solution of the set of differential equations
[t,y] = collocation(fname,0,1,y0,yf,guess,n,rho,[],w,fth);
% Temperature changes
for k = 1:n+2
options=optimset;
theta(k) = fzero(fth,30,options,y(:,k),w);
end
% Plotting the results
subplot(2,2,1), plot(t,y(1,:))
xlabel(‘Time’)
ylabel(‘Cell’)
title(‘(a)’)
subplot(2,2,2), plot(t,y(2,:))
xlabel(‘Time’)
ylabel(‘Penicillin’)
title(‘(b)’)
subplot(2,2,3), plot(t,y(3,:))
xlabel(‘Time’)
ylabel(‘First Adjoint’)
title(‘(c)’)
subplot(2,2,4), plot(t,theta)
xlabel(‘Time’)
ylabel(‘Temperature (deg C)’)
title(‘(d)’)

% Set dY for derivation
for k = 1:m*(n+1)
if Y(position(k)) ~= 0
dY(position(k)) = Y(position(k)) / 100;
else
dY(position(k)) = 0.01;
end
end

Cart

Cart