MATLAB: Using ode45 to solve a non-linear system of coupled ODE’s

coupledimplicitMATLABnonlinearodeode45system

so my code doesn't work. I've looked around and seen that maybe ode45 is not the best ode solver for this problem but it is an assignment and we're told to use that specific one.
So what we have is a system of 3 coupled first-order odes, where dx(1)/dt = vx = x(2), and dx(2)/dt = d(vx)/dt = ax = -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))
Very complicated function, unfortunately I'm struggling to simplify it but it shouldn't make too much of a difference hopefully!
So I'll post my code and the output, and hopefully someone can point me to a trivial mistake that I'm probably making. Keep in mind the problem I think lies with the fact that the coupled odes all depend on the other variables, so d(vx)/dt depends on vy = dy/dt and vz = dz/dt, for example.
Thanks in advance!
    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) [x(2); -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))];
odeset('abstol',1e-6);
[t,x] = ode45(odex,[0 time_length],[0,30/sqrt(2)]);
odey = @(t,y) [y(2); -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))];
odeset('abstol',1e-6);
[t,y] = ode45(odey,[0 time_length],[0,0]);
odez = @(t,z) [z(2); -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))];
odeset('abstol',1e-6);
[t,z] = ode45(odez,[0 time_length],[0,30/sqrt(2)]);
and error output:
Undefined function 'y' for input arguments of type 'double'.
Error in
@(t,x)[x(2);-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))]
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 24)
[t,x] = ode45(odex,[0 time_length],[0,30/sqrt(2)]);
I suspect I know a bit of what's going wrong, yet I'm at a loss for how to fix it. Any help is much appreciated!

Best Answer

  1. in your definition of W you're mixing doubles and function handles.
  2. write your ode-equations as one m-function - it becomes to messy with function handles.
  3. in that function you'll have something along:
function drdt = myode(t,r,whetever,more,you,need)
code
drdt(1) = vx;
drdr(2) = ax;
drdt(3) = vy;
drdr(4) = ay;
drdt(5) = vz;
drdr(6) = az;
then proceed as above.
Related Question