Need some help in solving these ODE's

Hey guys, I need to solve a coupled, simultaneous and implicit system of ODE's. Here's my code and then I'll briefly explain it:

close all
clear all

time_length = 25;
g = 9.8;             % Acceleration due to gravity (m/s^2)
rho0 = 1.20;         % Density of air (kg/m^3)
d = 0.22;            % Ball diameter (m)
m = 0.43;            % Ball mass (kg)
A = pi*(d/2)^2;      % Ball cross-sectional area (m^2)
Cd = 0.3;            % Drag coefficient
Cl = 0.3;            % Lift coefficient
Sx = 1;              % x-component of S
Sy = 1;              % y-component of S
Sz = 0;              % z-component of S
S = [Sx, Sy, Sz];    % Spin vector

% Coupled ODE's. x = x(1), vx = x(2). Same for y and z

      odex = @(t,x,y,z) -Cd*(A/(2*m))*rho0*sqrt(x(2).^2 + y(2).^2 + z(2).^2).*x(2) + Cl*(A/(2*m))*rho0*sqrt(x(2).^2 + y(2).^2 + z(2).^2).*(Sy.*z(2)-Sz.*y(2));

      odey = @(t,x,y,z) -Cd*(A/(2*m))*rho0*sqrt(x(2).^2 + y(2).^2 + z(2).^2).*y(2) + Cl*(A/(2*m))*rho0*sqrt(x(2).^2 + y(2).^2 + z(2).^2).*(Sz.*x(2)-Sx.*z(2));

      odez = @(t,x,y,z) -g-Cd*(A/(2*m))*rho0*sqrt(x(2).^2 + y(2).^2 + z(2).^2).*z(2) + Cl*(A/(2*m))*rho0*sqrt(x(2).^2 + y(2).^2 + z(2).^2).*(Sx.*y(2)-Sy.*x(2));

    W = @(x,y,z) [x(1);x(2);y(1);y(2);z(1);z(2)];

    dxyz = @(t,x,y,z) [x(2);odex;y(2);odey;z(2);odez];

    [T,W] = ode45(dxyz,linspace(0,time_length,1e4),[0,30/sqrt(2),0,0,30/sqrt(2)]);

So basically:

The notation is as follows:

x(1) = x,

x(2) = derivative of x = dx/dt = vx,

odex = second derivative of x = derivative of vx = d2x/dt2 = d(vx)/dt = dx(2)/dt = ax. Same for y and z.

Horrible equations in odex/y/z but basically each depend on the other variables, i.e. odex doesn't only depend on x, but y and z too (and their derivatives)

This last point means I need to solve them via matrix algebra. W is the column matrix that I'm planning to solve for, and dxyz is the column matrix containing the derivatives, i.e. dW/dt = dxyz.

in the ode45 bracket I have 6 initial conditions, corresponding to x(1) -> z(2) in order.

And now when I run it, I get these outputs:

Error using @(t,x,y,z)[x(2);odex;y(2);odey;z(2);odez]
Not enough input arguments.

Error in odearguments (line 88)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 114)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in code (line 37)
[T,W] = ode45(dxyz,linspace(0,time_length,1e4),[0,30/sqrt(2),0,0,30/sqrt(2)]);

It's telling me not enough input arguments, I do not understand. I have everything in there that I need, don't I?

Also I'm getting the idea that all these errors correspond to just the one error, i.e. if I fix the one thing that's wrong, I'll fix it all. Unfortunately, I don't know what to do about it because I don't understand why I'm getting an error.

Any and all help is much appreciated!

Answer 1

[Geoff Hayes]

Hi Joshua - while I haven't used ode45, I will try to give some hints at what may be causing the errors.

The message Not enough input arguments means that "someone" is not supplying enough input arguments to "some" function. In this case, that would be the ode45 not providing enough inputs to the dxyz function. Look at the documentation for ode45 and in particular the details on the function handle input parameter, odefun, which states

A function handle that evaluates the right side of the differential equations. All solvers solve systems of equations in the form y′ = f(t,y) or problems that involve a mass matrix, M(t,y)y′ = f(t,y). The ode23s solver can solve only equations with constant mass matrices. ode15s and ode23t can solve problems with a mass matrix that is singular, i.e., differential-algebraic equations (DAEs).

So the function must have a signature that expects two input parameters, t and y. Note how your function signature for dxyz requires four input parameters. Hence the error. So either your function signature will need to change or you will have to adapt the code to something similar to the second example (from the ode45 link).

Now note how dxyz is defined

 dxyz = @(t,x,y,z) [x(2);odex;y(2);odey;z(2);odez];

It produces a 6x1 vector, calling odex, odey, and odez. But the code is not supplying any input parameters to these functions so would generate an error. I observed the following errors

 Error using vertcat
 The following error occurred converting from double to function_handle:
 Error using function_handle
 Too many output arguments.

So you will need to change the above so that (whatever) inputs are passed into to these functions

 dxyz = @(t,x,y,z) [x(2);odex(t,x,y,z);y(2);odey(t,x,y,z);z(2);odez(t,x,y,z)];

As well, the initial condition vector, as supplied to ode45, has only five elements rather than six

 [0,30/sqrt(2),0,0,30/sqrt(2)]

-------------

The above addresses some concerns with the code, but not with whether it is being used as it should be. For example, why does the initial conditions vector have six elements when the three position elements (x,y,z) are unused? Why can't we just supply the velocities, and the dxyz return a 3x1 vector instead of the 6x1?

[Joshua D'Agostino]

Thanks for such a detailed answer, Geoff!

Okay, I completely get what you're saying about the input arguments, but the thing is that odex (say) is dependent on y and z too, making the inputs t, x, y and z, right? Is there a way around this?

Definitely have seen that vertcat error, I'll give odex/y/z input arguments then. That would require making a function for each wouldn't it? Because I actually tried that by doing something like this:

function odex(t,x,y,z)

     odex = { all that above }

end

function odey(t,x,y,z)
...

and so on, and I actually got error messages saying "A nested function cannot be defined inside a control statement" ... where is the control statement?

Good spot on the initial conditions, simple mistake!

[Geoff Hayes]

Joshua - you might be able to avoid the extra function handles (for odex, odey, and odez but creating a nested function for your dxyz. Usually, when I want to do something like this, I take my script (like your code above) and place it in a function where the name of the function is the name of your m-file

 function [T,W] = myFuncName

     close all;
     clear all;
     time_length = 25;

     % etc.

 end

The key thing to note is that all of your script code is "sandwiched" between function and end, with T and W being returned.

Now replace all of your odex, odey, odez, and dxyz code with the following nested function (note that I've kept the inputs as is until you decide how you want to change them

 function [out] = dxyz(t,v)

      x = v(1:2);
      y = v(3:4);
      z = v(5:6);

      odex = -Cd*(A/(2*m))*rho0*sqrt(x(2).^2 + y(2).^2 + z(2).^2).*x(2)...
         + Cl*(A/(2*m))*rho0*sqrt(x(2).^2 + y(2).^2 + z(2).^2).*(...
         Sy.*z(2)-Sz.*y(2));

      odey = -Cd*(A/(2*m))*rho0*sqrt(x(2).^2 + y(2).^2 + z(2).^2).*y(2)...
         + Cl*(A/(2*m))*rho0*sqrt(x(2).^2 + y(2).^2 + z(2).^2).*(...
         Sz.*x(2)-Sx.*z(2));

      odez =  -g-Cd*(A/(2*m))*rho0*sqrt(x(2).^2 + y(2).^2 + z(2).^2).*z(2)...
         + Cl*(A/(2*m))*rho0*sqrt(x(2).^2 + y(2).^2 + z(2).^2).*(...
         Sx.*y(2)-Sy.*x(2));

      out   = [x(2);odex;y(2);odey;z(2);odez];
 end

Following this, you can call the ode45 as before

 [T,W] = ode45(@dxyz,linspace(0,time_length,1e4),...
               [0,30/sqrt(2),0,0,30/sqrt(2) 42]);

Note the ampersand in front of the dxyz and the extra (dummy) element 42.

That should get the code to "work" without errors. The next step is to determine whether the solution makes sense, whether you really need six elements in your initial conditions and in the out from dxyz.

[Joshua D'Agostino]

Thanks heaps for your solution Geoff, and it does run without any issues. The only thing though is it seems to just be returning a time vector (1 to 25 in 10000 points), but nothing else. This isn't as expected, surely?

[Geoff Hayes]

T should be the time vector, but W should be the rest of the data. Make sure that when you call the function from the Command Window you do

 [T,W] = myFuncName;

When I ran the above, T was a 10000x1 vector and W was a 10000x6 vector.

[Joshua D'Agostino]

You sir, are a genius. Thank you so much for your help!

[Geoff Hayes]

Glad that I was able to help!

Answer 2

[Kate]

Hello, could I ask you relative questions here? I am just wondering if I can get the equations about the position in the x-,y-,z- direction? I have the data about the XYZposition, then I would like to get Cd and Cl from the positional data... Could you help me?

[Geoff Hayes]

MAMI - please post this as a question rather than as an answer.