# Need some help in solving these ODE's

19 views (last 30 days)
Joshua D'Agostino on 29 Sep 2014
Commented: Geoff Hayes on 12 Feb 2018
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!

Geoff Hayes on 29 Sep 2014
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?
Geoff Hayes on 1 Oct 2014
Glad that I was able to help!

Kate on 12 Feb 2018
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 on 12 Feb 2018
MAMI - please post this as a question rather than as an answer.