How to write the proper @ODEFUN for ode45?

G'day.

There is something funny with ode45 that if I change the arrangement of the element in the @ODE, the result change as well.

For example: Let write a function file:

function [f] = PosandVel(t,ic) %rho, Cd, R, and V are constant
%parameter
rho = 1000;
Cd = 1; %drag constant
M = 5;
R = 0.05;
V = 1; %velocity of fluid
fac = rho*Cd*pi*R^2/(2*M);

f(2)=ic(1); %dx/dt = u
f(1)=fac*(V^2-2*ic(1)+ic(1)^2); %du/dt=rho*Cd*pi*r/(2m)*(v^2-2uv+u^2)
f=f';

Then we write in the m-file:

clear
%Runge-Kutta method
t=[0:10]
ic=[0 0]
[T_OUT, ux] = ode45(@PosandVel,t,ic)

The result is as follows

ux =

         0         0
    0.4399    0.2620
    0.6110    0.7978
    0.7020    1.4583
    0.7585    2.1905
    0.7970    2.9694
    0.8249    3.7811
    0.8461    4.6170
    0.8627    5.4717
    0.8761    6.3413
    0.8870    7.2231

But if I change the element arrangement in the function file so that f(1) become f(2) and vice versa then the result become:

ux =

         0         0
         0    0.7854
         0    1.5708
         0    2.3562
         0    3.1416
         0    3.9270
         0    4.7124
         0    5.4978
         0    6.2832
         0    7.0686
         0    7.8540

Weird.

P.S. This is extracted from OCW MIT 2.29 Lecture 1.

http://ocw.mit.edu/courses/mechanical-engineering/2-29-numerical-fluid-mechanics-fall-2011/lecture-notes/

Answer 1

[Zhang lu]

the element arrangement cant be changed . ODE has strict rule.

f(1)=dx/dt=u, f(1) must be first-order differential equation

f(2)=dx2/dt2=du/dt=rho*Cd*pi*r/(2m)*(v^2-2uv+u^2), f(2) must be second-order differential equation

[Jan]

f(2) is not a 2nd order equation, but the 2nd component of the 1st order equation system. Matlab's ODE integrators require 1st order ODEs.

Answer 2

[Jan]

When you swap the components of the output of the ODE, you have to swap the inputs also. In your case this means, that ic(1) must be changed to ic(2) to obtain equivalent results.

With:

f(1) = ic(1);
f(2) = fac * (V ^ 2 - 2 * ic(1) + ic(1) ^ 2);

the first component cannot leave the position 0, because the derivative will be 0 also.

[Lydia]

Wow it does work! But since ic = [0 0], that is, both elements are zero anyway, why does changing the order have any effect?

[Jan]

@Lydia: f(1) is the derivative of the first component as f(2) belongs to the 2nd component. The derivative determines the slope of each component and is used to calculate the new values for the next time step. This new values are forwarded to the ODE as input "ic" in your case. So if "ic(1)" is zero and "f(1)=ic(1)", the first component start at zero and has the slope zero, such that it remains zero. But the second component of the derivative f(2) starts with "fac * (V ^ 2 -2*0 + 0^2)", which is not zero.

But when you swap the order of the derivatives, but not the order of the current positions "ic", you have another ODE. Not that the "ic" in the call of PosandVel equals the "ic" in the call of ODE45 only in the first step. Afterwards "ic" of PasandVel is updated.

I hope it is clearer now. If not, imagine that the effects are trivial. So drinking a coffee might reveal all details already.