How can I pass a large function into ODE45?

I have a large ODE function that I need to evaluate the system response time for. My function is a second order differential . In the code, theta_2 refers to the position of my link bar 2.

My problem is that my ODE is huge. Like really big. and

, --> , -->

, --> --> , -->

And is a peicewsie function of

What is the best way to pass my ODE into ODE45. Or any other ODE solver. My boundary conditions are

clc
clear
 

% Find time-response plots for theta2dot and theta 

% ----------- Knows -----------

R_2 = 2 ;           %inches
R_3 = 6 ;           %inches 
diam = 3.5;
radius = diam./2;

% ----------- Pressure Stuff ----------
n = 1.37 ;          % Gas R-value
P_atm = 14.7;
P1 = 13;            %psi 
P2 = P1;            %psi 
P3 = P2.*8.^n;      %psi 
P4 = 600;           %psi 
P5 = P4./(8.^n);    %psi 
P6 = 18;            %psi 
P7 = P6;            %psi 

% ------------- Input ------------- 

% theta_2 = 0:1:720;
% theta_2 = theta_2';

% Let theta_dot be omega 
% Let theta_2dot be alpha


R_4 = @(theta_2) R_2.*cosd(theta_2) + sqrt((R_3.^2)-(R_2*sind(theta_2)));

%R_4 = Get_R_4(theta_2);

theta_3 = @(theta_2) atan2d((-R_2).*sind(theta_2),(R_4-R_2.*cosd(theta_2)));

h_3 = @(theta_2) (-R_2.*cosd(theta_2))./(R_3.*cosd(theta_3));

h_3_prime =  @(theta_2)(R_2.*sind(theta_2))./(R_3.*cosd(theta_3)) + tand(theta_3).*(h_3.^2);

f_4 = @(theta_2) -R_2.*sind(theta_2) + R_2.*cosd(theta_2).*tand(theta_3);

f_4_prime = @(theta_2) -R_2.*cosd(theta_2) - R_3.*cosd(theta_3).*(h_3.^2) - R_3.*sind(theta_3).*h_3_prime;

R_4_prime = f_4;
R_4_2prime = f_4_prime;

% [pressure] = Get_pressure(theta_2);
R4_atBDC = R_4(180);
R4_atTDC = R_4(0);

d = 0.57142857;
        
BDC_Volume = (R4_atTDC-R4_atBDC+d).*pi.*radius.^2;
TDC_Volume = (d).*pi.*radius.^2;
Int_Volume = @(theta_2) ((R4_atTDC-R4_atBDC)-(R4(theta_2)-R4_atBDC)+d).*pi.*radius.^2;

Compression_Pressure = @(theta_2) P2.*(BDC_Volume).^1.37./(Int_Volume).^1.37;
Expansion_Pressure = @(theta_2) P4.*(TDC_Volume).^1.37./(Int_Volume).^1.37;

pressure = @(theta_2) piecewise( theta_2<180, P1, theta_2>180 & theta_2<360, ...
    Compression_Pressure, theta_2 == 360, P4, theta_2>360 & theta_2<540, ...
    Expansion_Pressure, theta_2>540 & theta_2<720, P6);


force = @(theta_2)(pressure(theta_2)-P_atm).*(radius.^2).*(pi);
Power = @(theta_2) f_4.*force;


Power_ODE = @(omega_2, alpha_2) Power == (0.02072+1.5+0.01166.*(R_4_prime.^2)).*(omega_2).*(alpha_2)+0.01166.*(R_4_prime).*(R_4_2prime).*(omega_2.^3)+0.005.*(omega_2.^3)+8.*cosd(theta_2).*(omega_2);


a = @(theta_2) (force(theta_2).*f_4(theta_2))./(0.02072+1.5+0.01.*R_4_prime(theta_2).^2);
b = @(theta_2) (-0.01166.*R_4_prime(theta_2).*R_4_2prime(theta_2))./(0.02072+1.5+0.01.*R_4_prime(theta_2).^2);
c = @(theta_2) (0.005)./(0.02072+1.5+0.01.*R_4_prime(theta_2).^2);
d = @(theta_2) 8./(0.02072+1.5+0.01.*R_4_prime(theta_2).^2);


odeFUNC = @(t,theta) [theta(2); a(theta(1)) - b(theta(1))*(theta(2)) - c(theta(1))*((theta(2))^2) - d(theta(1))*cosd(theta(1))];
tspan = [0 10];   
theta0 = [180 400];

[ts,ys] = ode45(@(t,theta) odeFUNC,tspan,theta0);

Error using superiorfloat
Inputs must be floats, namely single or double.

Error in odearguments (line 114)
dataType = superiorfloat(t0,y0,f0);

Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Answer 1

[Walter Roberson]

Your code needed a bunch of fix-ups. You defined a bunch of anonymous functions but then code them without passing any parameters to the anonymous functions. We can deduce the fix-ups for most (but not all!) of them from context.

You were using piecewise() on numeric values, but piecewise() is only defined for symbolic values.

% Find time-response plots for theta2dot and theta 

% ----------- Knows -----------

R_2 = 2 ;           %inches
R_3 = 6 ;           %inches R_4
diam = 3.5;
radius = diam./2;

% ----------- Pressure Stuff ----------
n = 1.37 ;          % Gas R-value
P_atm = 14.7;
P1 = 13;            %psi 
P2 = P1;            %psi 
P3 = P2.*8.^n;      %psi 
P4 = 600;           %psi 
P5 = P4./(8.^n);    %psi 
P6 = 18;            %psi 
P7 = P6;            %psi 

% ------------- Input ------------- 

% theta_2 = 0:1:720;
% theta_2 = theta_2';

% Let theta_dot be omega 
% Let theta_2dot be alpha


R_4_fh = @(theta_2) R_2.*cosd(theta_2) + sqrt((R_3.^2)-(R_2*sind(theta_2)));

%R_4_fh = Get_R_4(theta_2);

theta_3_fh = @(theta_2) atan2d((-R_2).*sind(theta_2),(R_4_fh(theta_2)-R_2.*cosd(theta_2)));

h_3_fh = @(theta_2) (-R_2.*cosd(theta_2))./(R_3.*cosd(theta_3_fh(theta_2)));

h_3_prime_fh =  @(theta_2)(R_2.*sind(theta_2))./(R_3.*cosd(theta_3_fh(theta_2))) + tand(theta_3_fh(theta_2)).*(h_3_fh(theta_2).^2);

f_4_fh = @(theta_2) -R_2.*sind(theta_2) + R_2.*cosd(theta_2).*tand(theta_3_fh(theta_2));

f_4_prime_fh = @(theta_2) -R_2.*cosd(theta_2) - R_3.*cosd(theta_3_fh(theta_2)).*(h_3_fh(theta_2).^2) - R_3.*sind(theta_3_fh(theta_2)).*h_3_prime_fh(theta_2);

R_4_prime_fh = f_4_fh;
R_4_2prime = f_4_prime_fh;

% [pressure_fh] = Get_pressure_fh(theta_2);
R4_atBDC = R_4_fh(180);
R4_atTDC = R_4_fh(0);

d_fh = 0.57142857;
        
BDC_Volume = (R4_atTDC-R4_atBDC+d_fh).*pi.*radius.^2;
TDC_Volume = (d_fh).*pi.*radius.^2;
INT_Volume_fh = @(theta_2) ((R4_atTDC-R4_atBDC)-(R_4_fh(theta_2)-R4_atBDC)+d_fh).*pi.*radius.^2;

Compression_Pressure = @(theta_2) P2.*(BDC_Volume).^1.37./(INT_Volume_fh(theta_2)).^1.37;
Expansion_Pressure = @(theta_2) P4.*(TDC_Volume).^1.37./(INT_Volume_fh(theta_2)).^1.37;

pressure_fh = @(theta_2) piecewise( theta_2<180, P1, theta_2>180 & theta_2<360, ...
    Compression_Pressure(theta_2), theta_2 == 360, P4, theta_2>360 & theta_2<540, ...
    Expansion_Pressure(theta_2), theta_2>540 & theta_2<720, P6);


force_fh = @(theta_2)(pressure_fh(theta_2)-P_atm).*(radius.^2).*(pi);
Power_fh = @(theta_2) f_4_fh(theta_2).*force_fh(theta_2);


Power_ODE = @(omega_2, alpha_2) Power_fh == (0.02072+1.5+0.01166.*(R_4_prime_fh.^2)).*(omega_2).*(alpha_2)+0.01166.*(R_4_prime_fh).*(R_4_2prime).*(omega_2.^3)+0.005.*(omega_2.^3)+8.*cosd(theta_2).*(omega_2);  %ERRROR Power_fh, R_4_prime_fh, R_4_2prime


a_fh = @(theta_2) (force_fh(theta_2).*f_4_fh(theta_2))./(0.02072+1.5+0.01.*R_4_prime_fh(theta_2).^2);
b_fh = @(theta_2) (-0.01166.*R_4_prime_fh(theta_2).*R_4_2prime(theta_2))./(0.02072+1.5+0.01.*R_4_prime_fh(theta_2).^2);
c_fh = @(theta_2) (0.005)./(0.02072+1.5+0.01.*R_4_prime_fh(theta_2).^2);
d_fh = @(theta_2) 8./(0.02072+1.5+0.01.*R_4_prime_fh(theta_2).^2);


odeFUNC = @(t,theta) [theta(2); a_fh(theta(1)) - b_fh(theta(1))*(theta(2)) - c_fh(theta(1))*((theta(2))^2) - d_fh(theta(1))*cosd(theta(1))];
tspan = [0 10];   
theta0 = [180 400];

syms T THETA [2 1]
FUNs = odeFUNC(T, THETA);
FUNh = matlabFunction(FUNs, 'vars', {T, THETA}, 'file', 'FUNh.m', 'optimize', false);

[ts,ys] = ode45(FUNh,tspan,theta0);
plot(ts, ys)
legend({'y1', 'y2'})

The plot appears to be empty because the outputs are constants, always [180, 400]

Your Power_ODE is wrong in three different places. You define it in terms of Power, R_4_prime and R_4_2prime, each of which are function handles defined in terms of theta_2, but you do not pass anything to any of them in Power_ODE, and we have no reason to guess whether omega_2 or alpha_2 should be what is passed. However it turns out not to matter because you do not invoke Power_ODE after defining it. You should probably also be defining Power_ODE in terms of piecewise()

Note that when you have something defined in terms of piecewise() then when you use matlabFunction() to generate an anonymous function, then you must use the 'file' output option.

Note that ode45() and all related ode*() functions are only valid if at all locations, the second derivatives of the expressions are continuous. In the great majority of the cases if you use piecewise expressions (or use interp1() with default linear interpolation) the first and second derivatives are not matched at the breakpoints, and you will either be given an error about failure to integrate, or else you will have the bad luck that you will not receive any error message. The reason not receiving any error message is bad luck is that you will not be notified that the answer you have gotten is garbage, and you are likely to trust it.

You do not happen to encounter any breakpoints in your particular ODE evaluated with those initial conditions and that time scale, so this particular problem did not immediately apply to you; as soon as you fix the expressions or initial conditions to do something more meaningful, then you will encounter problems.

I mentioned earlier that there were a lot of places that you define anonymous functions but then do not pass appropriate arguments to them. To help track that down, I renamed all of your anonymous functions to end with '_fh' to make it easier to read the expressions to make sure that the anonymous functions were being invoked with appropriate parameters each time.

I recommend that you adopt naming conventions that help you to reduce this kind of problem

[Walter Roberson]

Generally speaking, I advise you to mostly avoid mixing symbolic expressions (or symbolic functions) and anonymous functions. If you are needing to construct symbolically (perhaps you need constructs such as piecewise), then it is usually better to define everything symbolically until you get to the point of needing numeric evaluation, in which case you use matlabFunction() at that point.

Answer 2

[dpb]

Write your function as an m file and pass a handle to it; don't try to do stuff in scripts or complicated stuff as anonymous functions. Amongst other benefits, you can then much more easily test that it does as you expect at any point in time by calling it standalone with any input time and whatever auxiliary parameters are needed and debug it.

If you need additional parameters beyond the time and y' vector, see the section on <Parameterizing ODE Functions>

linked to in the doc for <ode45>