matlab invalid expression error

How do I correct this? I'm writing code for a 2nd order bvp, and have ann invalid expression error on this line.

function dydx = odefun( z, y, wc, (pi*D^2)/4*k, qs, {piDh)*(T-Ta) )
 
dydx = [ y(2);
         ( wc * y(2) - qs * y(1) - (piDh)*(T-Ta) ) / piD^2/4*k ];
 
end

I've tried adjusting the parentheses, not really sure what next to do.

[Walter Roberson]

When you use that syntax, what do you expect MATLAB to do with that 4th parameter, (pi*D^2)/4*k ? And is your expectation that MATLAB will consider that an exact match for piD^2/4*k even though the () are different ?

[izuchukwu morah]

the purpose is for matlab to map out a bvp using the necessary parameters involved in the base equation. (pi*D^2)/4*k is a parameter, more precisely, a coefficient from the equation.

and yes, the expectation would be for an exact match. Issue is, when I remove the parentheses, i get an error for improper use of "^", and it only goes away when i put back the parentheses which then switches to the first error i mentioned.

[Walter Roberson]

If I understand correctly, your expectation for

function dydx = odefun( z, y, wc, (pi*D^2)/4*k, qs, (piDh)*(T-Ta) )

is that anywhere in the code that (pi*D^2)/4*k occurs, you would want that expression replaced by the value passed in as the fourth parameter, and that anywhere in the code that (piDh)*(T-Ta) occurs, you would want that replaced by the value passed in as the sixth parameter?

What would your expectation be if the code contains

temp = (pi*D^2)/2*k

which, mathematically, would be twice the (pi*D^2)/4*k expression ? Would you expect it to be replaced with twice the value that was passed in as the fourth parameter? If the code contained (pi*D^2)/4 which is k times the value passed in as the fourth parameter, then what would you expect MATLAB to use as the value?

Your code does not currently contain (pi*D^2)/4*k and instead contains piD^2/4*k which is different not only in the () but in the fact that it contains a variable piD whereas the expression in the header was pi*D

If, hypothetically, you had had

function dydx = odefun( z, y, wc, (pi*D^2)/4*k, qs, (piDh)*(T-Ta), x^2+y^2 )

then is your expectation that MATLAB should be able to deduce the value of x^2-y^2 if that appeared in the code? Since, theoretically, it could see that it had x^2+y^2 and also had y, and so could square the y passed in and subtract twice that from the x^2+y^2 passed in in order to get the x^2-y^2 ? And if the expression in the code were x by itself, then hypothetically it could take the passed in x^2+y^2 and subtract the square of the passed in y, leaving x^2 which it could take the square root of... well, except for the fact that it would have lost the sign of x along the way...

[izuchukwu morah]

yes, those are my expectations. I still get the error even when i match up the ()

Answer 1

[Walter Roberson]

Please check that I extract the correct simulated output. You imply the first output should be extracted but the plot here is from the second output: the first output was around 4e5 for most of the range.

bvpmodel

piD24k = 628.3185

Elapsed time is 0.473840 seconds.

ans = 1×2

1 50

ans = 1×2

1 50

ans = 1×2

1 50

function bvpmodel

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                        %
% Inputs                                                                 %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% problem parameters: ====================================================
D=2;  % try: 1,3 
k=200; %try: 800
% coefficients:
wc  = 40.0;    % try: 30, 40, 50, 60
piD24k = pi*D^2/4*k
piD24k = 0.009;    % try: 0.006, 0.004, 0.0009    
qs = 0.0;    % try: 0, 1000, 10000
piDhTTa = 3.0;    % try: 32300, 6282, 8000, 10000    %piDhTTa is pi*Dh*(T-Ta)

% domain:
L  = 2;

% boundary conditions:
T0 = 30;
Ttarget = 200;

% solver parameters: =====================================================

% solution domain:
N  = 50;                 % number of nodes
z  = linspace( 0, L, N ); % set z equally spaced over the z domain

% initial guess: (constant values)
Tinit    = T0;            % T(z) = T0
dTdzinit =  0;            % dT/dz(z) = 0


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                        %
% Solution                                                               %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% solution initialization:
solinit = bvpinit( z, [ Tinit; dTdzinit ] );

% solution process - bvp4c solver:
% - the "@" notation is used in function handles; the expression:
%   "@( var1, var2, ... )fun( var1, var2, ..., vari, ... varn )"
%   allows using variables var1, var2, ... in calls to the function "fun" 
%   while the values of variables vari ... varn specified prior to the
%   call to bvp4c
%piD24k is pi*D^2/4*k
%piDhTTa is pi*Dh*(T-Ta)

tic  % start clock for computing time

sol = bvp4c( @( z, y )odefun( z, y, wc, piD24k, qs, piDhTTa),...
             @( ya, yb )bcfun( ya, yb, T0, Ttarget ),...
             solinit );
         
toc  % end clock for computing time


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                        %
% Post-Processing                                                        %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% retrieve solution:
%T = sol.y( 1, :);
temp = deval(sol, z);
T = temp(2,:);

% analytical solution for qs = 0:
Tanalytic = ( T0 + (Ttarget - T0) .* ( ( exp( wc.*z/(piD24k) ) - 1 )/( exp( wc.*L / (piD24k) ) - 1 ) ) + ( (piDhTTa) ).*L / (piD24k) ) .* ( ( z/L - exp( wc.*z / (piD24k) ) - 1 )/( exp( wc.*L / (piD24k) ) - 1 ) ) ;
%piD24k is pi*D^2/4*k
%piDhTTa is pi*Dh*(T-Ta)

% plots:
figure;
size(z), size(T), size(Tanalytic)
plot( z, T, '-b', z, Tanalytic, '+r');
grid;
box on;
legend( ' T ', ' T_{analytic} (for q_s = 0) ' );
title(' Boundary Value Problem ' );
xlabel( ' z ' );
ylabel( ' T(z) ');

end

function dydx = odefun( z, y, wc, piD24k, qs, piDhTTA )
%piD24k is pi*D^2/4*k
%piDhTTa is pi*Dh*(T-Ta)
dydx = [ y(2);
         ( wc * y(2) - qs * y(1) - piDhTTA ) / piD24k ];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                        %
% Function evaluating the boundary conditions: res = F( ya, yb, ... )    %
%                                                                        %
% For this problem: y = [ y1 y2 ], y1 = T, y2 = dT/dz                    %
%                                                                        %
% x = 0: T - T0 = 0 --> y1 - T0 = 0                                      %
% x = L: T - Ttarget = 0 --> y1 - Ttarget = 0                                      %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function res = bcfun( ya, yb, T0, Ttarget )
res = [ ya(1) - T0  ;
        yb(1) - Ttarget ];
    
end

[izuchukwu morah]

i undertsand what you did, the simulated output should plot from 0.1 on z-axis, up till T0=30 on T(z) axis. similar to this sample.

Answer 2

[Walter Roberson]

MATLAB does not work according to your expectations. No computer languages reason about the relationships between variables automatically to deduce values not stated. With the possible exception of SNOBOL or similar languages.

function dydx = odefun( z, y, wc, piD24k, qs, piDhTTa )
%piD24k is pi*D^2/4*k
%piDhTTa is pi*Dh*(T-Ta)
 
dydx = [ y(2);
         ( wc * y(2) - qs * y(1) - piDhTTa ) / piD24k ];
 
end

See, give full names for the expressions and use the same names where you want the values to be substituted. No deduction needs to be done to find the values of the variables. Use better variable names if you want.

[izuchukwu morah]

i understand what you mean, that helps a little, but random errors keep showing up. now it's error for variable Tanalytic. (full code)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
% Solution of a Boundary Value Problem (BVP) using Matlab.
%
% The problem is given by:
% wc * dT/dz - (piD^2/4*k) * d^2T/dz^2 + qs - (pi*Dh)(T-Ta)=0,   0 > z > L
%
% The boundary conditions for the problem are:
%
%             @ z = 0, T = T0
%             @ z = L: T = Ttarget
%
% For qs = 0, the problem has the analytical solution:
%
%  T(z) = T0 + (Ttarget - T0) * ( ( exp( wc*z / (pi * d^2/4k) ) - 1 )/( exp( wc*L / (pi*d^2/4k) ) - 1 ) ) ...
%   + ( (pi*D*h)(T-Ta) )*L / (pi*D^2 / 4k) ) * ( ( z/L - exp( wc*z / (pi*D^2 / 4k) ) - 1 )/( exp( wc*L / (pi*D^2 / 4k) ) - 1 ) )
%
% Solution of BVP using Matlab: ==========================================
%
% Solution of the boundary value problem: dy/dx = odefun( z, y )
% with boundary conditions:   res( za, ya ) = 0, res( zb, yb ) = 0
% and using as initial guess: yinit
% 
% % solution initialization:
% solinit = bvpinit( x, yinit );
% 
% % solver options:
% options = bvpset( 'RelTol', 1e-5 );
%
% % solver call:
% sol = bvp4c( @odefun, @bcfun, solinit, options );
% 
% % required functions:
% dydx = odefun( z, y );
% res  = bcfun( ya, yb );
% 
% % output:
% sol.z     nodes of the mesh selected by bvp4c
% sol.y      approximation to y(z) at the mesh points of sol.z
% sol.yp      approximation to dy/dz(z) at the mesh points of sol.z
% sol.solver  type of solver used, i.e. 'bvp4c'
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function bvpmodel

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                        %
% Inputs                                                                 %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% problem parameters: ====================================================

% coefficients:
wc  = 40.0;    % try: 30, 40, 50, 60
piD24k  = 0.009;    % try: 0.006, 0.004, 0.0009    %piD24k is pi*D^2/4*k
qs = 0.0;    % try: 0, 1000, 10000
piDhTTa = 3.0;    % try: 32300, 6282, 8000, 10000    %piDhTTa is pi*Dh*(T-Ta)

% domain:
L  = 2;

% boundary conditions:
T0 = 30;
Ttarget = 200;

% solver parameters: =====================================================

% solution domain:
N  = 100;                 % number of nodes
x  = linspace( 0, L, N ); % set x equally spaced over the x domain

% initial guess: (constant values)
Tinit    = T0;            % T(z) = T0
dTdzinit =  0;            % dT/dz(z) = 0


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                        %
% Solution                                                               %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% solution initialization:
solinit = bvpinit( x, [ Tinit; dTdzinit ] );

% solution process - bvp4c solver:
% - the "@" notation is used in function handles; the expression:
%   "@( var1, var2, ... )fun( var1, var2, ..., vari, ... varn )"
%   allows using variables var1, var2, ... in calls to the function "fun" 
%   while the values of variables vari ... varn specified prior to the
%   call to bvp4c
%piD24k is pi*D^2/4*k
%piDhTTa is pi*Dh*(T-Ta)

tic  % start clock for computing time

sol = bvp4c( @( z, y )odefun( z, y, wc, piD24k, qs, piDhTTa),...
             @( ya, yb )bcfun( ya, yb, T0, Ttarget ),...
             solinit );
         
toc  % end clock for computing time


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                        %
% Post-Processing                                                        %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% retrieve solution:
T = sol.y( 1, : );

% analytical solution for qs = 0:
% Tanalytic = T0 + (Ttarget - T0) * ( ( exp( wc*z / (pid24k) ) - 1 )/( exp( wc*L / (pid24k) ) - 1 ) ) ...
%   + ( (piDhTTa) )*L / (piD24k) ) * ( ( z/L - exp( wc*z / (piD24k) ) - 1 )/( exp( wc*L / (piD24k) ) - 1 ) ) ;
%piD24k is pi*D^2/4*k
%piDhTTa is pi*Dh*(T-Ta)
% plots:
figure;
plot( x, T, '-b', x, Tanalytic, '+r');
grid;
box on;
legend( ' T ', ' T_{analytic} (for q_s = 0) ' );
title(' Boundary Value Problem ' );
xlabel( ' z ' );
ylabel( ' T(z) ');

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                        %
% Function describing the BVP: dy/dx = F( x, y )                         %
%                                                                        %
% For this problem: y = [ y1; y2 ], y1 = T, y2 = dT/dz                   %
%                                                                        %
%          wc * dT/dz - (piD^2/4*k) * d^2T/dz^2 + qs - (pi*Dh)(T-Ta)=0            %
%                                                                        % 
% => F  = [ F1; F2 ]:                                                    %
%    F1 = dy1/dz = y2 = dT/dz                                            %
%    F2 = dy2/dz = d^2T/dz^2 = ( wc * y2 - qs * y1 - (piDh)(T-Ta) )/(piD^2/4k)        %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function dydx = odefun( z, y, wc, piD24k, qs, piDhTTA )

%piD24k is pi*D^2/4*k
%piDhTTa is pi*Dh*(T-Ta)

dydx = [ y(2);
         ( wc * y(2) - qs * y(1) - piDhTTA ) / piD24k ];

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                        %
% Function evaluating the boundary conditions: res = F( ya, yb, ... )    %
%                                                                        %
% For this problem: y = [ y1 y2 ], y1 = T, y2 = dT/dz                    %
%                                                                        %
% x = 0: T - T0 = 0 --> y1 - T0 = 0                                      %
% x = L: T - Ttarget = 0 --> y1 - Ttarget = 0                                      %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function res = bcfun( ya, yb, T0, Ttarget )

res = [ ya(1) - T0  ;
        yb(1) - Ttarget ];
    
end

[Steven Lord]

These are the only places in your code where Tanalytic appears.

% Tanalytic = T0 + (Ttarget - T0) * ( ( exp( wc*z / (pid24k) ) - 1 )/( exp( wc*L / (pid24k) ) - 1 ) ) ...
%   + ( (piDhTTa) )*L / (piD24k) ) * ( ( z/L - exp( wc*z / (piD24k) ) - 1 )/( exp( wc*L / (piD24k) ) - 1 ) ) ;
%piD24k is pi*D^2/4*k
%piDhTTa is pi*Dh*(T-Ta)
% plots:
figure;
plot( x, T, '-b', x, Tanalytic, '+r');

You never compute the value of the variable Tanalytic (the line of code that would do so is commented out) so that variable doesn't exist on the line where you try to plot it.

[Walter Roberson]

Your code never assigns to that variable, except in some comments.

[izuchukwu morah]

okay, i see what both of you mean! last error that pops up though is a vector error when i try to plot. I'm not sure where to match up the vectors for 'T' and 'Tanalytic'

function bvpmodel

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                        %
% Inputs                                                                 %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% problem parameters: ====================================================
D=2;  % try: 1,3 
k=200; %try: 800
% coefficients:
wc  = 40.0;    % try: 30, 40, 50, 60
piD24k = pi*D^2/4*k
piD24k = 0.009;    % try: 0.006, 0.004, 0.0009    
qs = 0.0;    % try: 0, 1000, 10000
piDhTTa = 3.0;    % try: 32300, 6282, 8000, 10000    %piDhTTa is pi*Dh*(T-Ta)

% domain:
L  = 2;

% boundary conditions:
T0 = 30;
Ttarget = 200;

% solver parameters: =====================================================

% solution domain:
N  = 50;                 % number of nodes
z  = linspace( 0, L, N ); % set z equally spaced over the z domain

% initial guess: (constant values)
Tinit    = T0;            % T(z) = T0
dTdzinit =  0;            % dT/dz(z) = 0


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                        %
% Solution                                                               %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% solution initialization:
solinit = bvpinit( z, [ Tinit; dTdzinit ] );

% solution process - bvp4c solver:
% - the "@" notation is used in function handles; the expression:
%   "@( var1, var2, ... )fun( var1, var2, ..., vari, ... varn )"
%   allows using variables var1, var2, ... in calls to the function "fun" 
%   while the values of variables vari ... varn specified prior to the
%   call to bvp4c
%piD24k is pi*D^2/4*k
%piDhTTa is pi*Dh*(T-Ta)

tic  % start clock for computing time

sol = bvp4c( @( z, y )odefun( z, y, wc, piD24k, qs, piDhTTa),...
             @( ya, yb )bcfun( ya, yb, T0, Ttarget ),...
             solinit );
         
toc  % end clock for computing time


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                        %
% Post-Processing                                                        %
%                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% retrieve solution:
T = sol.y( 1, :);

% analytical solution for qs = 0:
Tanalytic = ( T0 + (Ttarget - T0) .* ( ( exp( wc.*z/(piD24k) ) - 1 )/( exp( wc.*L / (piD24k) ) - 1 ) ) + ( (piDhTTa) ).*L / (piD24k) ) .* ( ( z/L - exp( wc.*z / (piD24k) ) - 1 )/( exp( wc.*L / (piD24k) ) - 1 ) ) ;
%piD24k is pi*D^2/4*k
%piDhTTa is pi*Dh*(T-Ta)

% plots:
figure;
plot( z, T, '-b', z, Tanalytic, '+r');
grid;
box on;
legend( ' T ', ' T_{analytic} (for q_s = 0) ' );
title(' Boundary Value Problem ' );
xlabel( ' z ' );
ylabel( ' T(z) ');

end

Answer 3

[Steven Lord]

When you define your function, the function declaration line should include the names of the variables into which the input arguments will be stored. It cannot contain an expression that is not a name.

When you call your function, specify the exact values on which you want your function to operate.

This is an incorrect way to define an addme function:

function z = addme(2, 3)
z = 2+3;
end

This is a correct way to define the addme function:

function z = addme(x, y)
z = x + y;
end

and this is a correct way to call the addme function.

theOutput = addme(2, 3)

See this documentation page for more information on defining functions.