Solving Coupled Second Order ODE by ode45

Hi, I'm trying to solve the 2nd order ODE's

S1'' = -k * sqrt( S1'^2 + S2'^2 ) * S1',
S2'' = -k * sqrt( S1'^2 + S2'^2 ) * S2' - g.

Which I have attempted to do by writing it as four first order ODE's

V1' = -k * sqrt( V1^2 + V2^2 ) * V1,
V2' = -k * sqrt( V1^2 + V2^2 ) * V2 - g,
S1' = V1,
S2' = V2.

To solve this I've tried to use ODE45. Originally I solved just V1 and V2 which worked fine, by making a function

function vp = Func(t,v)

k = 1
g = 2;
vp = diag( [ - k .* sqrt( v .^ 2 + v(2) .^ 2 ) .* v, - k .* sqrt( v(1) .^ 2 + v .^ 2 ) .* v - g] );
end

t0 = 0;
tfinal = 120;
v0 = [ 1 2 ]'
tfinal = tfinal*(1+eps);
tspan = t0:dt:tfinal;
[t,v] = ode45(@Func,tspan,v0);

which works fine. I thought this would be a simple extension to solving the four differentials at once by changing this to

function vp = Func(t,v)

k = 1
g = 2;
vp = diag( [ - k .* sqrt( v .^ 2 + v(2) .^ 2 ) .* v, - k .* sqrt( v(1) .^ 2 + v .^ 2 ) .* v - g, v(1), v(2)] );
end

t0 = 0;
tfinal = 120;
v0 = [ 1 2 3 4]'
tfinal = tfinal*(1+eps);
tspan = t0:dt:tfinal;
[t,v] = ode45(@Func,tspan,v0);

However, when I try this I get the error Dimensions of matrices being concatenated are not consistent. Any hints in the right direction as to how to extend from solving a system of two differential equations to four would be greatly appreciated.

[Star Strider]

I’m not quite sure what you’re doing. You’ve created ‘vp’ as a (6x6) diagonal matrix. It might work better as a column vector.

[McTavish Dylan]

Hi, thanks very much for the response. I created a diagonal matrix as to write my ODE in the form M(t,y)y′ = f(t,y) as the documentation for ode45 says to do.

I'm a bit confused that you say it's a 6x6 matrix though, surely it should be 4x4?

Thanks again

[Star Strider]

I deleted the code, but when I substituted ‘v’ as an arbitrary (1x2) vector, it created a (6x6) diagonal matrix for ‘vp’.

Actually, the MATLAB ODE functions (at least all of them I’ve worked with thus far) want the derivatives to be a column vector.

Meanwhile, see if the code in my Answer helps.

[McTavish Dylan]

I can't try your code right at the moment (don't have access to a computer right now) but I will as soon as I do.

Thanks very much once again for the prompt help.

[McTavish Dylan]

Hi, I don't quite understand how I would adapt your answer that uses syms to use ode45, however using your advice to write it as a column vector I have managed to get it to work correctly, so thank you for solving my problem.

[Star Strider]

My pleasure!

My Answer gives the full Symbolic Math Toolbox derivation. You would use the ‘Sys’ anonymous function for your ODE function, with appropriate changes to get it to work with your code.

So it would magickally transform to your ‘Func’ function as:

Sys = @(Y,g,k)[Y(2);-g-k.*sqrt(Y(2).^2+Y(4).^2).*Y(2);Y(4);-k.*sqrt(Y(2).^2+Y(4).^2).*Y(4)];
Func = @(t,y) Sys(Y,g,k);

with ‘g’ and ‘k’ previously defined in your workspace.

[Ryan Compton]

Hi, I am trying to implement this code as a skeleton to adapt to my own data. However, running it as is (changing the vector from diagonal to column) I get the following error, I am not sure how this is even possible:

Error using odearguments (line 95) FUNC returns a vector of length 10, but the length of initial conditions vector is 4. The vector returned by FUNC and the initial conditions vector must have the same number of elements.

Did you run into this issue? I know you said you were able to get it to work.

Answer 1

[Star Strider]

Not yet being awake enough to do this on paper, I let the Symbolic Math Toolbox loose with your system:

syms S1(t) S2(t) k g
% S1'' = -k * sqrt( S1'^2 + S2'^2 ) * S1',
% S2'' = -k * sqrt( S1'^2 + S2'^2 ) * S2' - g.
EqS1 = diff(S1,2) == -k * sqrt( diff(S1)^2 + diff(S2)^2 ) * diff(S1);
EqS2 = diff(S2,2) == -k * sqrt( diff(S1)^2 + diff(S2)^2 ) * diff(S2) - g;
[ODE,Vars] = odeToVectorField(EqS1, EqS2)
Sys = matlabFunction(ODE)

ODE =
                                 Y[2]
 - g - k*(Y[2]^2 + Y[4]^2)^(1/2)*Y[2]
                                 Y[4]
      -k*(Y[2]^2 + Y[4]^2)^(1/2)*Y[4]

Vars =
  S2
 DS2
  S1
 DS1

Sys = @(Y,g,k)[Y(2);-g-k.*sqrt(Y(2).^2+Y(4).^2).*Y(2);Y(4);-k.*sqrt(Y(2).^2+Y(4).^2).*Y(4)]

Make appropriate changes in ‘Sys’ to create a function that ode45 will like. See if that does what you want.