Solving Coupled Second Order ODE by ode45

3 次查看（过去 30 天）
McTavish Dylan
McTavish Dylan 2015-12-20
评论： Ryan Compton 2018-9-25
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.
  7 个评论
显示 5更早的评论
Star Strider
Star Strider 2015-12-20
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
Ryan Compton 2018-9-25
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.

请先登录，再进行评论。

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

采纳的回答

Star Strider
Star Strider 2015-12-20
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.

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by