How do I pass inputs into external functions from ode45?

I'm trying to run an ode45 function, but the differential equations are constrained by other equations, which I made functions for, and referenced in the ode function (func) , but I keep getting the error "Not enough input arguments." There may be other errors, but I can't get past this one. I need the domain of the ode45 (pi to 3pi) to be passed into the other functions as input theta. My main goal is to plot theta against pressure/work but I can't get it to intialize. Thanks in advance

global r gamma thetas thetab l S b V0 T1 Tw Qin R 
gamma = 1.4;
r = 8.4;
thetas = (3*pi)/2;
thetab = pi;

% meters
l = 0.012;
S = 0.08;
b = 0.09;

% cm squared
V0 = 50;

% Kelvin 
T1 = 300;
Tw = 300;

% J/kg
Qin = 2.8e6;

% J/kg/k
R = 287;

% sol = ode45(@func,[pi 3*pi],[1.013e5 0]);
[theta,sol] = ode45(@(theta,y) func,[pi 3*pi],[1.013e5 0]);

function referenced by ode45 and all the dependent functions

function E2 = func(theta,y)
global gamma b Tw Qin R
%y(1)=Pressure
%y(2)=Work

%mass blow by
he = 0;

%heat transfer to cylinder wall
% C = 0;

%instantaneous surface area
Aw =(4*V(theta))/b;

%temperature
T = (y(1)*V(theta))/(mass(theta)*R);

%Pressure and work derivatives, respectively
E2 = [-gamma*(y(1)/V(theta))*(V_prime(theta))+(gamma-1)*((mass(theta)*Qin)/V(theta))*(x_prime(theta))-((gamma-1)*he*Aw*(T-Tw))/(V*w)+((gamma*mass_prime(theta))/mass(theta))*y(1);
        y(1)*V_prime(theta);];

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Volume function
function Vtheta = V(theta)
global r l S
sigma = S/(2*l);
Vtheta = 1 + ( (r-1) / (2*sigma) )*( 1+sigma*(1-cos(theta)-sqrt(1-sigma^2*sin(theta)^2)));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% x derivative function
function dxdtheta = x_prime(theta)

global thetas thetab

if pi <= theta && theta < thetas
    dxdtheta = 0;
elseif thetas <= theta && theta <= (thetas+thetab)
    dxdtheta = -(pi/( 2* thetab))*cos( ( pi*( theta-thetas ) ) / thetab);

elseif thetas+thetab < theta && theta < 3*pi
    dxdtheta = 0;

end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%mass function
function M = mass(theta)

M0 = 1;%initial mass
C =0 ;%heat transfer to cylinder
omega = 1;%rotational velocity
M = M0*exp( (-C/omega) * (theta-pi) );
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Mass derivative
function dMdTheta = mass_prime(theta)

%constants
omega_prime = 0;%rotational acceleration - zero bc/ angular velocity was zeero
M0 = 1;%initial mass
C =0 ;%heat transfer to cylinder

dMdTheta = ( -(C*omega_prime)-(pi*omega_prime) )*M0*exp( (-C/omega)*(theta-pi) );

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Volume derivative
function dVdTheta = V_prime(theta)
global r l S
sigma = S/(2*l);

dVdTheta = ( 2*r*sqrt( 1-( sigma^2 ) * sin( theta )^2 )* sin( theta ) + r*sigma*sin(2*theta)-2*sqrt( 1-(sigma^2)*( sin(theta) )^2 )*sin( theta ) )/( 4*sqrt( 1 - (sigma^2)*( sin( theta ) )^2 ) );

end

[Walter Roberson]

See

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

Note that

[theta,sol] = ode45(@(theta,y) func,[pi 3*pi],[1.013e5 0]);

tells ode45 that it will be operating on a function that accepts two parameters, and ignores them, and returns whatever func() with no parameter would return.

We recommend that you get rid of the global variables and parameterize the functions instead.

[Joshua]

Okay so besides parametrizing, would writing it like this be better?

sol = ode45(@func,[pi 3*pi],[1.013e5 0]);

[Walter Roberson]

If you are not going to parameterize to get rid of the globals, then Yes, @func would be faster than @(theta,y)func(theta,y)

[Joshua]

I'd like to parametrize them but im not sure how to in this context. Would you take a look at the code in my response to Torstens answer and give me some feedback?

Answer 1

[Walter Roberson]

Parameterized version -- no globals.

Note: I had to guess what V*w meant.

p.gamma = 1.4;
p.r = 8.4;
p.thetas = (3*pi)/2;
p.thetab = pi;

% meters
p.l = 0.012;
p.S = 0.08;
p.b = 0.09;

% cm squared
p.V0 = 50;

% Kelvin 
p.T1 = 300;
p.Tw = 300;

% J/kg
p.Qin = 2.8e6;

% J/kg/k
p.R = 287;

%other
p.omega = 1;%rotational velocity
p.omega_prime = 0;%rotational acceleration - zero bc/ angular velocity was zeero

% sol = ode45(@func,[pi 3*pi],[1.013e5 0]);
[theta,sol] = ode45(@(theta,y) func(theta,y,p),[pi 3*pi],[1.013e5 0]);

Warning: Failure at t=6.040709e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.421085e-14) at time t.

function E2 = func(theta,y,p)
%y(1)=Pressure
%y(2)=Work

%mass blow by
he = 0;

%heat transfer to cylinder wall
% C = 0;

%instantaneous surface area
Aw =(4*V(theta,p))/p.b;

%temperature
T = (y(1)*V(theta,p))/(mass(theta,p)*p.R);

%Pressure and work derivatives, respectively
E2 = [-p.gamma*(y(1)/V(theta,p))*(V_prime(theta,p))+(p.gamma-1)*((mass(theta,p)*p.Qin)/V(theta,p))*(x_prime(theta,p))-((p.gamma-1)*he*Aw*(T-p.Tw))/(V(theta,p)*p.Tw)+((p.gamma*mass_prime(theta,p))/mass(theta,p))*y(1);
        y(1)*V_prime(theta,p);];

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Volume function
function Vtheta = V(theta,p)
sigma = p.S/(2*p.l);
Vtheta = 1 + ( (p.r-1) / (2*sigma) )*( 1+sigma*(1-cos(theta)-sqrt(1-sigma^2*sin(theta)^2)));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% x derivative function
function dxdtheta = x_prime(theta,p)
if pi <= theta && theta < p.thetas
    dxdtheta = 0;
elseif p.thetas <= theta && theta <= (p.thetas+p.thetab)
    dxdtheta = -(pi/( 2* p.thetab))*cos( ( pi*( theta-p.thetas ) ) / p.thetab);
elseif p.thetas+p.thetab < theta && theta < 3*pi
    dxdtheta = 0;
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%mass function
function M = mass(theta,p)
M0 = 1;%initial mass
C =0 ;%heat transfer to cylinder
M = M0*exp( (-C/p.omega) * (theta-pi) );
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Mass derivative
function dMdTheta = mass_prime(theta,p)

%constants


M0 = 1;%initial mass
C =0 ;%heat transfer to cylinder

dMdTheta = ( -(C*p.omega_prime)-(pi*p.omega_prime) )*M0*exp( (-C/p.omega)*(theta-pi) );

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Volume derivative
function dVdTheta = V_prime(theta,p)
sigma = p.S/(2*p.l);

dVdTheta = ( 2*p.r*sqrt( 1-( sigma^2 ) * sin( theta )^2 )* sin( theta ) + p.r*sigma*sin(2*theta)-2*sqrt( 1-(sigma^2)*( sin(theta) )^2 )*sin( theta ) )/( 4*sqrt( 1 - (sigma^2)*( sin( theta ) )^2 ) );

end

[Joshua]

Thanks so much, the V*w was an error that I should have already fixed. Is there any way to get rid of the warning about integration tolerances, or is that entirely dependent on the equations Im using?

[Walter Roberson]

I had to guess about what V*w was intended to mean, and I might have guessed wrong, so I might have made the equations mostly unsolvable over the range you want.

But if I did happen to guess correctly about what it was intended to mean, then you have a problem that your current equations lead to a singularity at just over 6 seconds.

if pi <= theta && theta < p.thetas

The mathematics of ode45 and related functions requires that the second derivatives of the functions be continuous during any one call to ode45. Any time you see an if statement in your code that calculates the function, you should strongly suspect that you have failed to provide continuous second derivatives. Your code should probably be creating event functions for the two differerent crossings, with the event functions configured to terminate the function, after which you would restart with the new conditions.

Answer 2

[Torsten]

This code works without syntax errors, but you should pass all your parameters to "func" and distribute them from there among the other functions called instead of using globals.

Take a look here on how to do this:

https://de.mathworks.com/help/matlab/math/parameterizing-functions.html

global r gamma thetas thetab l S b V0 T1 Tw Qin R 
gamma = 1.4;
r = 8.4;
thetas = (3*pi)/2;
thetab = pi;

% meters
l = 0.012;
S = 0.08;
b = 0.09;

% cm squared
V0 = 50;

% Kelvin 
T1 = 300;
Tw = 300;

% J/kg
Qin = 2.8e6;

% J/kg/k
R = 287;

% sol = ode45(@func,[pi 3*pi],[1.013e5 0]);
[theta,sol] = ode45(@func,[pi 3*pi],[1.013e5 0]);

Warning: Failure at t=6.040709e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.421085e-14) at time t.

plot(theta,sol)

Warning: Imaginary parts of complex X and/or Y arguments ignored.

function referenced by ode45 and all the dependent functions

function E2 = func(theta,y)
global gamma b Tw Qin R
%y(1)=Pressure
%y(2)=Work

%mass blow by
he = 0;

%heat transfer to cylinder wall
% C = 0;

%instantaneous surface area
Aw =(4*V(theta))/b;

%temperature
T = (y(1)*V(theta))/(mass(theta)*R);

%Pressure and work derivatives, respectively
E2 = [-gamma*(y(1)/V(theta))*(V_prime(theta))+(gamma-1)*((mass(theta)*Qin)/V(theta))*(x_prime(theta))-((gamma-1)*he*Aw*(T-Tw))/(V(theta)*Aw)+((gamma*mass_prime(theta))/mass(theta))*y(1);
        y(1)*V_prime(theta);];

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Volume function
function Vtheta = V(theta)
global r l S
sigma = S/(2*l);
Vtheta = 1 + ( (r-1) / (2*sigma) )*( 1+sigma*(1-cos(theta)-sqrt(1-sigma^2*sin(theta)^2)));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% x derivative function
function dxdtheta = x_prime(theta)

global thetas thetab

if pi <= theta && theta < thetas
    dxdtheta = 0;
elseif thetas <= theta && theta <= (thetas+thetab)
    dxdtheta = -(pi/( 2* thetab))*cos( ( pi*( theta-thetas ) ) / thetab);

elseif thetas+thetab < theta && theta < 3*pi
    dxdtheta = 0;

end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%mass function
function M = mass(theta)

M0 = 1;%initial mass
C =0 ;%heat transfer to cylinder
omega = 1;%rotational velocity
M = M0*exp( (-C/omega) * (theta-pi) );
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Mass derivative
function dMdTheta = mass_prime(theta)

%constants
omega_prime = 0;%rotational acceleration - zero bc/ angular velocity was zeero
omega=1;
M0 = 1;%initial mass
C =0 ;%heat transfer to cylinder

dMdTheta = ( -(C*omega_prime)-(pi*omega_prime) )*M0*exp( (-C/omega)*(theta-pi) );

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Volume derivative
function dVdTheta = V_prime(theta)
global r l S
sigma = S/(2*l);

dVdTheta = ( 2*r*sqrt( 1-( sigma^2 ) * sin( theta )^2 )* sin( theta ) + r*sigma*sin(2*theta)-2*sqrt( 1-(sigma^2)*( sin(theta) )^2 )*sin( theta ) )/( 4*sqrt( 1 - (sigma^2)*( sin( theta ) )^2 ) );

end

[Joshua]

I see that you changed the ode45 arguments that I had, but did you make other changes? I got rid of the global variables and defined the constants inside the "func" function, as well as nesting the other functions in "func" but I'm still getting the "Not enough input arguments." error. I have a feeling Im not correctly parametrizing this but I'm not sure how to proceed further.

% sol = ode45(@func,[pi 3*pi],[1.013e5 0]);
[theta,sol] = ode45(@func,[pi 3*pi],[1.013e5 0]);


%New attempt at making an equation to be used with ODE45
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function E2 = func(theta,y)
gamma = 1.4;
r = 8.4;
thetas = (3*pi)/2;
thetab = pi;

% meters
l = 0.012;
S = 0.08;
b = 0.09;

% cm squared
V0 = 50;

% Kelvin
T1 = 300;
Tw = 300;

% J/kg
Qin = 2.8e6;

% J/kg/k
R = 287;
%y(1)=Pressure
%y(2)=Work

%mass blow by
he = 0;

%heat transfer to cylinder wall
% C = 0;

%instantaneous surface area
Aw =(4*V(theta))/b;

%temperature
T = (y(1)*V(theta))/(mass(theta)*R);

%Pressure and work derivatives, respectively
E2 = [-gamma*(y(1)/V(theta))*(V_prime(theta))+(gamma-1)*((mass(theta)*Qin)/V(theta))*(x_prime(theta))-((gamma-1)*he*Aw*(T-Tw))/(V*w)+((gamma*mass_prime(theta))/mass(theta))*y(1);
    y(1)*V_prime(theta);];


%Volume function
    function Vtheta = V(theta)
       
        sigma = S/(2*l);
        Vtheta = 1 + ( (r-1) / (2*sigma) )*( 1+sigma*(1-cos(theta)-sqrt(1-sigma^2*sin(theta)^2)));
    end

% x derivative function
    function dxdtheta = x_prime(theta)

        if pi <= theta && theta < thetas
            dxdtheta = 0;
        elseif thetas <= theta && theta <= (thetas+thetab)
            dxdtheta = -(pi/( 2* thetab))*cos( ( pi*( theta-thetas ) ) / thetab);

        elseif thetas+thetab < theta && theta < 3*pi
            dxdtheta = 0;

        end
    end

%mass function
    function M = mass(theta)

        M0 = 1;%initial mass
        C =0 ;%heat transfer to cylinder
        omega = 1;%rotational velocity
        M = M0*exp( (-C/omega) * (theta-pi) );
    end

%Mass derivative
    function dMdTheta = mass_prime(theta)

        %constants
        omega_prime = 0;%rotational acceleration - zero bc/ angular velocity was zeero
        M0 = 1;%initial mass
        C =0 ;%heat transfer to cylinder

        dMdTheta = ( -(C*omega_prime)-(pi*omega_prime) )*M0*exp( (-C/omega)*(theta-pi) );

    end

%Volume derivative
    function dVdTheta = V_prime(theta)
      
        sigma = S/(2*l);

        dVdTheta = ( 2*r*sqrt( 1-( sigma^2 ) * sin( theta )^2 )* sin( theta ) + r*sigma*sin(2*theta)-2*sqrt( 1-(sigma^2)*( sin(theta) )^2 )*sin( theta ) )/( 4*sqrt( 1 - (sigma^2)*( sin( theta ) )^2 ) );

    end

end

[Torsten]

Both V and w in the product V*w are undefined in E2.

omega is undefined in mass_prime.

[Walter Roberson]

You have

E2 = [-gamma*(y(1)/V(theta))*(V_prime(theta))+(gamma-1)*((mass(theta)*Qin)/V(theta))*(x_prime(theta))-((gamma-1)*he*Aw*(T-Tw))/(V*w)+((gamma*mass_prime(theta))/mass(theta))*y(1);
    y(1)*V_prime(theta);];

Notice the

/(V*w)

That is a problem in that statement because V is not a variable: it is a function of theta.

Also, there is no w anywhere in your code -- only Tw

[Joshua]

I see. It seems that was an issue from the beginning. Probably should have checked for that before posting this. Thank you for your help, I really appreciate it.