using ode45 to solve a non-linear system of coupled ODE's

Question

Joshua D'Agostino 2014-9-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/156525-using-ode45-to-solve-a-non-linear-system-of-coupled-ode-s

回答： Bjorn Gustavsson ，2014-10-5

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!

1 个评论
显示无隐藏无

Joshua D'Agostino 2014-9-28

bump

请先登录，再进行评论。

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

Answer 1

Bjorn Gustavsson 2014-10-5

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/156525-using-ode45-to-solve-a-non-linear-system-of-coupled-ode-s#answer_154013

in your definition of W you're mixing doubles and function handles.
write your ode-equations as one m-function - it becomes to messy with function handles.
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.

0 个评论
显示 -1更早的评论隐藏 -1更早的评论

请先登录，再进行评论。

Answer 2

Jan 2014-9-28

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/156525-using-ode45-to-solve-a-non-linear-system-of-coupled-ode-s#answer_153181

The error message seems to be clear:

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))]

This is a function of the arguments t and x, but the contents of the function uses the variables x, y and z.

What are y and z?

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

Joshua D'Agostino 2014-9-29

yeah right I recognised that, and so previously I had given them all @(t,x,y,z) variable handles. This doesn't solve anything though. New error message:

Error using
@(t,x,y,z)[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))]
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 24)
[t,x] = ode45(odex,[0 time_length],[0,30/sqrt(2)]);

Do you understand what I'm trying to do though? odex for example is the coupled ode's in the x-direction, where odex(1) = dx/dt and odex(2) = d(vx)/dt = d2x/dt2. It just so happens that all 3 odex, odey and odez are dependent on all x, y and z. That's why I initially just had @(t,x) on odex, because with that one I am only concerned with the x-direction. It's clear why it doesn't work, but I don't know, is there a way I can solve all 3 simultaneously using ode45? Because correct me if I'm wrong, but that would solve the problem, wouldn't it?

Joshua D'Agostino 2014-9-30

okay cheers I've gone and done that, but what is this giving me now and why?

    odex = @(t,x) -Cd*(A/(2*m))*rho0*sqrt(x(2).^2 + x(4).^2 + x(6).^2).*x(2) + Cl*(A/(2*m))*rho0*sqrt(x(2).^2 + x(4).^2 + x(6).^2).*(Sy.*x(6)-Sz.*x(4));
    odey = @(t,x) -Cd*(A/(2*m))*rho0*sqrt(x(2).^2 + x(4).^2 + x(6).^2).*x(4) + Cl*(A/(2*m))*rho0*sqrt(x(2).^2 + x(4).^2 + x(6).^2).*(Sz.*x(2)-Sx.*x(6));
    odez = @(t,x) -g-Cd*(A/(2*m))*rho0*sqrt(x(2).^2 + x(4).^2 + x(6).^2).*x(6) + Cl*(A/(2*m))*rho0*sqrt(x(2).^2 + x(4).^2 + x(6).^2).*(Sx.*x(4)-Sy.*x(2));
W = @(t,x) [x(1);x(2);x(3);x(4);x(5);x(6)];
dxyz = @(t,x) [x(2);@out;x(4);odey;x(6);odez];
[T,W] = ode45(dxyz,linspace(0,time_length,1e4),[0,30/sqrt(2),0,0,0,30/sqrt(2)]);

output:

Error using vertcat
The following error occurred converting from double to function_handle:
Error using function_handle
Too many output arguments.
Error in @(t,x)[x(2);@out;x(4);odey;x(6);odez]
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 35)
[T,W] = ode45(dxyz,linspace(0,time_length,1e4),[0,30/sqrt(2),0,0,0,30/sqrt(2)]);

I don't quite understand why it's giving me these errors, I thought I've taken care of the number of input arguments?

Jan 2014-10-5

I do not understand your code. What is "@out"? Why do you define W as a function, but overwrite it be the output of the integrator? "odey" is a function handle. If you concatenate it with "[x(2); odey]" what do you expect as result? Which type should this vector have?

I suggest to create a simple function instead of a pile of nested function handles, which are obviously too confusing.

请先登录，再进行评论。