How do I pass inputs into external functions from ode45?

Question

0 个投票

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

Walter Roberson 2023-4-28

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 2023-4-28

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?

请先登录，再进行评论。

请先登录，再回答此问题。

Follow Question

Answer 1

Walter Roberson 2023-4-29

在 MATLAB Online 中打开

0 个投票

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

2 个评论
显示无隐藏无

Joshua 2023-4-29

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 2023-4-29

在 MATLAB Online 中打开

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 2023-4-28

编辑：Torsten 2023-4-28

在 MATLAB Online 中打开

0 个投票

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

Joshua 2023-4-28

在 MATLAB Online 中打开

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

Walter Roberson 2023-4-29

编辑：Walter Roberson 2023-4-29

在 MATLAB Online 中打开

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 2023-4-29

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.

请先登录，再进行评论。

How do I pass inputs into external functions from ode45?

4 个评论
显示 2更早的评论隐藏 2更早的评论

采纳的回答

2 个评论
显示无隐藏无

更多回答（1 个）

4 个评论
显示 2更早的评论隐藏 2更早的评论

类别

产品

版本

标签

Community Treasure Hunt

How do I pass inputs into external functions from ode45?

4 个评论 显示 2更早的评论 隐藏 2更早的评论

采纳的回答

2 个评论 显示 无 隐藏 无

更多回答（1 个）

4 个评论 显示 2更早的评论 隐藏 2更早的评论

类别

产品

版本

标签

另请参阅

Community Treasure Hunt

4 个评论
显示 2更早的评论隐藏 2更早的评论

2 个评论
显示无隐藏无

4 个评论
显示 2更早的评论隐藏 2更早的评论