What is wrong with my code, my equation is vdv/dy=-G*M/(y+R)^2 ?

3 次查看(过去 30 天)
cssCopy codefunction dydt = gravity(y, v, G, M, R)
% y: height above surface
% v: velocity at height y
% G: gravitational constant
% M: mass of the planet
% R: radius of the planet
dydt(1,1) = v; % v = dy/dt
dydt(2,1) = -G*M/(y+R)^2; % a = d^2y/dt^2 = dv/dt
end
scssCopy codeG = 6.67430e-11; % gravitational constant
M = 5.9722e24; % mass of Earth
R = 6.371e6; % radius of Earth
  7 个评论
Dyuman Joshi
Dyuman Joshi 2023-4-4
编辑:Dyuman Joshi 2023-4-4
@Laura I am well aware that you are trying to solve the equation by MATLAB, that is not my point.
For what time range are you trying to solve the equation? What are the initial conditions to the differential equation?

请先登录,再进行评论。

回答(2 个)

Sam Chak
Sam Chak 2023-4-4
编辑:Sam Chak 2023-4-4
This is a simple example of using the solver's syntax that you can find in the ode45() documentation. If there are universal constants that they don't change, then place them inside the ode function, as suggested by @Dyuman Joshi. Otherwise, parameters that you would like to test with different value (i.e., mass), then you can pass the parameters to the ode function like @Cris LaPierre did.
What you need to do is to replace my free fall ODE with your description of Newton's law of universal gravitation.
Logically, when the object (point mass) hits the ground, then the simulation should stop. Anyhow, get familiar with the basic code. Click into the documentation for more examples.
% solving the ODE
tspan = [0 6]; % simulation time
y0 = [100 0]; % initial condition y(1) = 100 m; y(2) = 0 m/s
[t, y] = ode45(@freefall, tspan, y0);
% plotting the solution y(1)
plot(t, y(:,1), 'linewidth', 1.5),
xlabel('Time (sec)'), ylabel('Height (m)')
yline(0, '--', 'Ground level');
yregion(-80, 0)
% Describing the ODE
function dydt = freefall(t, y)
m = 1; % point mass
g = 9.8203; % gravity
dydt(1,1) = y(2); % dy/dt
dydt(2,1) = -g/m; % assume no air drag
end

Torsten
Torsten 2023-4-4
syms G M R y v(y)
Dv = diff(v,y);
eqn = diff(v,y,2) == -G*M/(y+R)^2;
cond = [v(0)==100,Dv(0)==0];
sol = dsolve(eqn,cond)
sol = 
Gnum = 6.67430e-11; % gravitational constant
Mnum = 5.9722e24; % mass of Earth
Rnum = 6.371e6; % radius of Earth
sol = vpa(subs(sol,[G M R],[Gnum,Mnum,Rnum]))
sol = 
fplot(sol,[0 6])
  2 个评论
Dyuman Joshi
Dyuman Joshi 2023-4-4
@Torsten, the equation you wrote is incorrect.
It's
v*diff(v,y)
instead of
diff(v,y,2)
And it needs only 1 intial condition
Torsten
Torsten 2023-4-4
编辑:Torsten 2023-4-4
I also thought so first , but Laura's code (and all answers and comments thereafter) speak another language.
I think d^2v/dy^2 was misinterpreted as d(0.5*v^2)/dy.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Programming 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by