フィルターのクリア

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

2 ビュー (過去 30 日間)
Lily
Lily 2024 年 4 月 23 日
編集済み: Sam Chak 2024 年 4 月 23 日
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!!

採用された回答

Steven Lord
Steven Lord 2024 年 4 月 23 日
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
  1 件のコメント
Lily
Lily 2024 年 4 月 23 日
Okay, I think I understand. Thank you so much for the help!

サインインしてコメントする。

その他の回答 (2 件)

Star Strider
Star Strider 2024 年 4 月 23 日
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')
.
  4 件のコメント
Steven Lord
Steven Lord 2024 年 4 月 23 日
I'd leave this for others (who are allowed to use Symbolic Math Toolbox) to learn from.
Star Strider
Star Strider 2024 年 4 月 23 日
@Steven Lord — O.K.

サインインしてコメントする。


Sam Chak
Sam Chak 2024 年 4 月 23 日
編集済み: Sam Chak 2024 年 4 月 23 日
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')

カテゴリ

Help Center および File ExchangeOrdinary Differential Equations についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by