How can I write a set of first order ODEs with two variables?

Hi, my assignment asks me to convert the following two 2nd order ODE's to a system of four 1st order ODEs and then encode the equations in a function to then use with ODE45.

The two 2nd order ODE's are as follows:

Then I converted them to a system of four 1st order ODEs (I think I did this correctly):

y' = v

x' = u

u' = -(D/m)*u*sqrt(u^2+v^2)

v' = -g-(D/m)*v*sqrt(u^2+v^2)

Then, I've been looking at this page as a reference for writing the system in my code. But, I'm confused as to how to write it when there are two different equations that use both x and y. The problem does say to slit the equations into several first order ODE's in terms of the values of y. Thank you!!

Answer 1

[Steven Lord]

In this case, the second input to your ODE function will have four elements. Using a slight variant (to avoid using y in two different ways) of the notation from that documentation page, vec′=f(t,vec), we have that vec is [x; y; u; v]. [The order doesn't really matter; you could have vec be [y; x; u; v] as well if that's more convenient. As long as you keep the same order throughout your entire code it will work.] That means you need your function f to return [x'; y'; u'; v'] (which is vec').

What is x'? From your transformation, that's just u which is the third element of the vec input to your function f.

What is y'? Again from the transformation, that's just v which is the fourth element of the vec input to your function f.

Similarly, you can compute u' and v' using the elements from vec as you've done mathematically.

Here's the general outline, you just need to fill in the FILL_THIS_IN sections. Then use this function when you call the ODE solver.

function vecprime = f(t, vec)
% Unpack vec into its components, so the expressions for the *prime variables
% will "look like" the mathematical form of the ODE system
x = vec(1);
y = vec(2);
u = vec(3);
v = vec(4);
% Define the components of vecprime
xprime = u;
yprime = v;
uprime = FILL_THIS_IN;
vprime = FILL_THIS_IN;
% Pack the *prime variables into vecprime
vecprime = [xprime; yprime; uprime; vprime];    
end

[Lily]

Okay, I think I understand. Thank you so much for the help!

Answer 2

[Star Strider]

Try something like this —

syms D g m t x(t) y(t) Y  
dx = diff(x);
d2x = diff(dx);
dy = diff(y);
d2y = diff(dy);

DEqn1 = d2x == -(D/m)*dx*sqrt(dx^2 + dy^2)

DEqn1(t) = 

DEqn2 = d2y == -g -(D/m)*dy*sqrt(dx^2 + dy^2)

DEqn2(t) = 

[VF,Subs] = odeToVectorField(DEqn1, DEqn2)

VF = 

Subs = 

DE12fcn = matlabFunction(VF, 'Vars',{t,Y,D,g,m})

DE12fcn = function_handle with value:
    @(t,Y,D,g,m)[Y(2);-g-(D.*sqrt(Y(2).^2+Y(4).^2).*Y(2))./m;Y(4);-(D.*sqrt(Y(2).^2+Y(4).^2).*Y(4))./m]

D = rand

D = 0.4206

m = 42;
g = 9.81;

ic = [1 0.1 1 0.1];

[t,y] = ode45(@(t,y)DE12fcn(t,y,D,g,m), [0 1], ic);

figure
plot(t, y)
grid
xlabel('t')
ylabel('Value')
legend(string(Subs), 'Location','best')

.

[Lily]

Thank you for the response. I'm not sure how much this impacts your answer, but I'm not allowed to use syms.

[Star Strider]

That negates it.

I’ll delete it in a few minutes.

[Steven Lord]

I'd leave this for others (who are allowed to use Symbolic Math Toolbox) to learn from.

[Star Strider]

@Steven Lord — O.K.

Answer 3

[Sam Chak]

Hi @Lily

I made some edits to enhance the annotations. Essentially, the idea is to represent the dynamic system in state-space format. Here is another approach you can try the ode45 solver. This approach also provides the flexibility to use the latest SUNDIALS Solvers, such as cvodesNonstiff and cvodesStiff.

%% Setup ODE problem
F               = ode;              % create an 'ode' object
F.Parameters    = [1 2 3];          % D = 1, m = 2, g = 3
F.ODEFcn        = @(t, x, p) [x(3);                                             % x' = u
                              x(4);                                             % y' = v
                           - (p(1)/p(2))*x(3)*sqrt(x(3)^2 + x(4)^2);            % u' = ...
                           - (p(1)/p(2))*x(4)*sqrt(x(3)^2 + x(4)^2) - p(3)];    % v' = ...
F.InitialValue  = [9; 7; 5; 3];     % x(0) = 9, y(0) = 7, u(0) = 5, v(0) = 3
F.Solver        = "ode45";          % can also select the lastest "SUNDIALS" Solvers

%% Ready to solve the ODE
tStart          = 0;
tFinal          = 10;
sol             = solve(F, tStart, tFinal);

%% Plot results
plot(sol.Time, sol.Solution,"-o"), grid on, xlabel('t')
title('System responses')
legend('x(t)', 'y(t)', 'u(t)', 'v(t)', 'location', 'southwest')